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TASK  OBJECTIVES 


There  are  two  broad  objectives  in  this  research  effort.  The  first  objective  is  to  improve 
our  understanding  of  the  propagation  of  high  frequency  regional  phases  through  the 
continental  crust.  This  analysis  involves  the  synthesis  of  realistic  wavetrains,  array 
processing  and  data  analysis  of  NORESS/PINESS,  Soviet  and  U.S.  data.  The  second 
objective  is  to  investigate  the  generation  of  seismic  noise  at  the  seafloor  and, 
subsequently,  the  continents. 

Specifically,  the  Work  Statement  includes  the  following; 


•  Compare  regional  broadband  seismic  data  from  the  following  sources: 

-  NORESS/PINESS 

-  Soviet  high  frequency  data  from  the  IRIS/IDA  network  in  the  USSR. 

-  U.S.  broad  band  stations  including  the  Pihon  IRIS/IDA  station  and  the 
stations  established  in  the  western  US  for  monitoring  tests  in  Nevada  and 
monitored  at  IGPP  through  a  satellite  telemetry  link. 

•  Synthesize  regional  phases  through  the  use  of  wavenumber  integration. 

-  Examine  complex  crustal  and  upper  mantle  structures. 

-  Examine  Q  arxl  frequency  dependence  at  the  European,  Soviet  and  US 
sites. 

•  Array  Processing  of  NORESS/PINESS  data  and  comparison  to  synthetics. 

-  Determine  the  frequency  dependence  of  polarization  anomalies  caused  by 
anisotropy.  Apply  to  NORESS,  Soviet  and  U.S.  data. 

-  Investigate  the  role  of  scattering  in  generation  of  regional  phase  coda. 

•  Physical  modeling  of  seafloor  and  continental  noise. 
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TECHNICAL  RESULTS 


ABSTRACT 

In  this  report  we  present  the  results  of  investigations  dealing  with  three  topics: 

•  We  have  continued  our  investigation  of  the  influence  the  regular  repetition  of 
features  in  seisniograms,  introduced  at  the  source  and/or  during  propagation  by 
layer  resonance,  has  on  the  spectrum  of  the  recorded  coda.  We  have  developed  an 
automatic  algorithm  that  discriminates  mining  explosions  from  nuclear 
explosions  and  Earthquakes. 

•  We  consider  the  influence  large  crustal  inhomogeneities,  or  topographic 
undulations,  in  the  vicinity  of  receivers  can  have  on  the  character  of  the  recorded 
seismic  coda.  In  specific  we  have  developed  a  migration  technique  that  scans  the 
recorded  coda  for  phases  generated  by  local  scattering  interactions. 

•  We  have  developed  a  regularized  approach  to  seismic  deconvolution  that  allows 
the  introduction  of  a  priori  information,  such  as  that  describing  the  statistical 
nature  of  the  additive  noise  and  the  underlying  model,  into  the  inversion  process. 
We  investigate  how  this  algorithm  can  enhance  the  images  produced  by  the 
migration  algorithm  developed  in  the  previous  section. 

SUMMARY 

As  part  of  the  last  DARPA/AFGL  contract  we  endeavored  to  determine  how  regional 
seismic  network  data  can  be  used  to  discriminate  ripple  fired  mining  explosions  from 
single  event  explosions.  As  part  of  this  work,  theory  was  developed  which  makes  two 
predictions.  Ripple-fired  quarry  blasts  should  produce  coda  that  is  highly  colored.  This 
spectral  character  should  be  resistant  to  change  with  time.  A  technique  (sonogram 
analysis)  was  developed  which  expands  the  information  contained  in  time  series  into  a 
time-frequency  plane  and  thus  allows  visual  identification  of  time-independant  spectral 
components.  The  technique  was  originally  applied  to  data  collected  by  the  NRDC  network 
of  sensors.  This  dataset  proved  to  be  well  suited  to  this  endeavor  since  it  contains 
regional  recordings  of  single-event  chemical  explosions  (detonated  in  the  fall  of  1987 
by  a  joint  Soviet-American  team)  in  addition  to  a  large  number  of  recordings  of  events 
strongly  believed  to  be  ripple-fired  mining  explosions.  The  proximity  of  the  stations  to 
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the  events  is  important  since  the  algorithm  depends  on  (and  exploits)  the  retention  of 
high-frequency  energy  -  it  is  a  broad  band  technique. 

This  work  has  been  extended  in  Chapter  1  to  use  data  recorded  by  arrays  (in  this  case  the 
NOR  ESS  small  aperture  array  in  Norway)  and  to  determine  if  ripple-fired  events  in 
other  geographic,  and  technical,  settings  appear  to  be  anomalous  when  viewed  using  the 
sonogram  technique.  In  view  of  the  inverse  relationship  between  the  size  a  seismic  event 
and  its  frequency  of  occurrence,  this  chapter  is  devoted  to  the  development  of  a  technique 
that  automatically  discriminates  between  quarry  blasts  and  earthquakes,  or  single-event 
explosions.  It  is  based  on  the  assumption  that  quarry  blasts  that  are  comprised  of  a 
number  of  sub  explosions  closely  grouped  in  space  and  time  can  be  discriminated  from 
the  simpler  events  by  virtue  of  the  fact  that  they  produce  coda  with  distinctive  spectral 
features  that  are  resistant  to  change  as  time  passes. 

In  chapter  2,  work  that  has  bearing  on  nuclear  detection  and  discrimination  of  relatively 
large  events  (roughly  150  kT)  from  teleseismic  ranges  is  introduced.  This  chapter 
describes  how  coda  recorded  by  arrays  can  be  scanned,  or  migrated,  to  image  the  local 
crust  for  scatterers  large  enough  to  emit  recognizable  surface  waves  (Rg)  when  excited 
by  teleseismic  body  phases.  There  is  some  concern  that  such  arrivals  may  be  confused 
with  the  pP  phase  used  to  discriminate  deep  events  from  shallow  events  (most 
earthquakes  from  nuclear  explosions).  Although  more  physically  comprehensive 
methods  (such  as  f-k  migration,  Kirchoff  migration  and  finite  difference  migration) 
have  been  developed  to  process  seismic  reflection  data,  we  have  chosen  to  adapt 
hyperbola-summation  migration  to  our  problem  because  we  are  dealing  with  dispersive 
wavelets.  The  idea  of  using  small-aperture  arrays  for  this  type  of  imaging  is  tested 
extensively  by  examining  synthetic  (and  thus  controlled)  datasets.  The  synthetic 
sources  are  infinitely  small  and  can  be  used  to  assess  the  spatial  resolution  of  the 
technique,  and  how  this  varies  as  the  source  moves,  the  array  changes,  or  the  imaging 
parameters  vary.  More  significantly,  this  chapter  concludes  by  relating  the  radial 
resolution  with  the  temporal  duration  of  the  teleseismic  wavelrains.  It  is  concluded  that 
deconvolution  of  the  array  records  is  necessary,  and  can  be  used  to  increase  radial 
resolution.  Unlike  the  work  in  Chapter  1.  which  requires  the  retention  of  high 
frequency  energy  and  thus  uses  regional  seismic  stations,  this  chapter  deals  with  a 
relatively  narrow,  and  low,  band  of  seismic  energy  -  from  1 .0  to  3.0  Hz. 


In  Chapter  3  a  technique  that  relies  on  inverse  theory  is  developed  to  deconvolve  the 
teleseismic  wavetrain  from  the  seismic  recordings  to  improve  the  radial  spatial 
resolution.  The  algorithm  is  designed  to  produce  models  by  trading  off  two  undesirable 
quantities  -  data  misfit  and  model  complexity.  A  priori  information,  obtained  from 
physical  arguments  regarding  how  we  expect  incident  seismic  energy  interacts  with 
large  crustal  scatterers,  is  used  to  define  what  level  of  model  simplicity,  in  the  form  of 
serial  correlation  or  smoothness,  should  be  present  in  the  model.  The  statistical 
character  of  the  additive  noise  (which  is  assumed  to  be  due  to  a  random  process)  is 
estimated  by  assuming  ergodicity  and  examining  pre-event  noise. 
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Chapter  1 


An  automatic  means  to  discriminate  between 
earthquakes  and  quarry  blasts 

Abstract 

In  this  chapter  we  discuss  our  efforts  to  use  the  NORESS  array  to  discriminate  between  regional 
earthquakes  and  ripple-fircd  quarry  blasts  (events  that  involve  a  number  of  sub-explosions  closely  grouped 
in  space  and  time).  The  method  we  describe  is  an  extension  of  the  lime  versus  frequency  “pattern-based" 
discriminant  proposed  by  Hedlin  et  al.  (1989b).  At  the  heart  of  the  discriminant  Ls  the  observation  that 
ripple-hred  events  tend  to  give  rise  to  coda  dominated  by  prominent  spectral  features  that  arc  independent 
of  time  and  periodic  in  frequency.  This  spectral  character  is  generally  absent  from  the  coda  produced 
by  earthquakes  and  “single-event"  explosions.  The  discriminant  originally  proposed  by  Hedlin  et  al 
(1989b)  used  data  collected  at  2.‘i0  s  '  by  single  sensors  in  the  1987  NRDC  network  in  Kazakhstan. 
U.S.S.R..  We  have  found  that  despite  the  relatively  low  digitization  rate  provided  by  the  NORESS  array 
(40  .s  ’)  we  have  had  good  success  in  our  efforts  to  discriminate  between  earthquakes  and  quarry  blasts 
by  stacking  all  vertical  array  channels  to  improve  signal  to  noise  ratios. 

We  describe  our  efforts  to  automate  the  method,  so  that  visual  pattern  recognition  is  not  required, 
and  to  make  it  less  susceptible  to  spurious  time-independent  spectral  features  not  originating  at  the 
source.  In  essence  we  compute  a  Fourier  transform  of  the  time-frequency  matrix  and  examine  the 
power  levels  representing  energy  that  is  periodic  in  frequency  and  independent  of  time.  Since  a  double 
Fourier  transform  is  involved,  our  method  can  be  considered  as  an  extension  of  "cepstral"  aruilysi<, 
(Tribotet,  1979).  We  have  found,  however,  that  our  approach  is  superior  since  it  is  cognizant  of  the  time 
independence  of  the  spcxiral  features  of  interest.  Wc  use  earthquakes  to  define  what  cepstral  power  is  in 
be  expected  in  the  absence  of  ripple  firing  and  search  for  events  that  violate  this  limit.  The  assessment 
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of  the  likelihood  that  ripple  tiring  occuned  at  the  source  is  made  automatically  by  the  computer  and  is 
based  on  the  extent  to  which  the  limit  is  violated. 

1.1  Introduction 

There  is  a  peculiar  breed  of  seismic  event  known  as  a  ripple-fired  explosion.  Such  an  event 
differs  markedly  from  a  standard  “single-event”  explosion  since  it  involves  tlie  detonation  of  numerous 
sub-explosions  closely,  and  generally  regularly,  grouped  in  space  and  time.  Ripple-tiring  is  a  technique 
commonly  u.sed  in  quarry  blasting  (Langefors  and  Kihlstrdm,  1978)  where  mine  operators  are  striving  to 
reduce  ground  motions  in  areas  proximal  to  the  mine,  enhance  rtxk-fracturing  and  reduce  the  amount  of 
material  thrown  into  the  air  -  “tly”  or  "throw”  rock  -  (Dowding,  1985).  Ripple-tiring  is  in  widespread 
u.sc,  being  employed  both  in  the  Americas  and  in  Europe  (Stump  et  al.,  1989). 

There  has  been  increased  interest  in  recent  years  in  discriminating  mining  events  from  earth¬ 
quakes  and  nuclear  explosions.  A  reduced  Threshold  Test  Ban  Treaty  could  potentially  bring  the  magni¬ 
tude  of  the  largest  nuclear  explivsions  down  to  that  of  large  "engineering"  explosions  otherwise  known  as 
quarry  blasts  (Stump  and  Reamer,  1988).  Aggravating  the  problem  is  the  existence  of  numerous  quarries 
in  the  vicinity  of  the  Scmipalatinsk  nuclear  lest  site  in  the  Soviet  Union  (Thurber  et  al.,  1989;  Hedlin 
t  t  111  I989h).  There  have  been  a  number  of  studies  dealing  dira'tly  and  indirectly  with  this  problem. 
Uxiking  primarily  at  Scandinavian  events  recorded  by  the  NORESS  array,  Baumgardi  and  Ziegler  (1988) 
luund  piominent  spectral  modulaliiHi  in  events  believed  to  involve  ripple-tiring,  but  not  in  the  spectra 
computed  from  earthquake  seismograms.  Hedlin  t’(  al  (1989b)  observed  similar  spectral  modulation  in 
tlic  cixla  produced  by  suspected  quarry  blasts  in  Ka/akhsi;m.  U.S.S.R.,  but  not  in  the  coda  produced  by 
single-event  calibration  explosions  detonated  at  similar  ranges.  They  found  further  that  the  modulation, 
when  present,  was  independent  of  time  from  the  onset,  well  into  llic  L,  cixla.  This  time-independeni 
character  has  also  been  observed  in  the  ccxla  prixluccd  by  quarry  blasts  and  recorded  in  Scandinavia 
(Hedlin  ei  al .  1989a).  Both  Baumgardt  and  Ziegler  (1988)  and  Hedlin  ei  al  (1989b)  found  that  the 
spectral  nuxlulatinn  observed  in  the  cixLi  produced  by  mine  explosions  could  be  rcprixluccd  effectively 
by  assuming  that  all  sub-expiosions  prixlucc  the  same,  common,  wavcfixm  and  that  the  motions  su¬ 
perpose  linearly.  Stump  and  Reinke  (1988)  have  investigated  the  validity  of  the  assumption  of  linear 
superposition.  They  prixluccd  strong  evidence  suppirrUng  the  assumption  when  wavclields  from  small, 
cliiscly  spaced,  explosions  are  ob.scrved  in  the  ncarficld.  Baumgardi  and  Ziegler  (1988).  Hedlin  et  al 
( I'lK'i.i.bi,  Slump  aixl  Reamer  (1988)  and  Smith  (1989)  -  who  also  observed  prominent  peaks  in  the 
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spectra  of  phases  produced  by  some  quarry  blasts  -  all  concluded  that  the  unusual  spectral  color  could 
be  used  to  discriminate  quarry  blasts  from  other  events  with  “whiter"  spectra. 

In  this  work  we  arc  extending  the  study  described  in  Hedlin  et  a/.(  1989b),  hereafter  referred  to 
as  paper  1,  in  a  number  of  ways.  We  examine  recordings  of  earthquakes,  not  single-event  explosions, 
to  determirte  if  they  can  be  discriminated  from  quarry  blasts  with  a  similar  degree  of  success.  We  seek 
to  determine  the  sensitivity  of  the  method  to  the  recording  “environment".  The  recordings  examined  in 
the  current  study  have  been  made  at  40  by  the  NORESS  small  aperlure-array  in  Norway  (Ringdal 
and  Husebye,  1982;  Mykkeltvcit  el  al.,  1983).  The  data  considered  in  paper  1  were  recorded  by  single 
sensors  and  digitized  at  250  s  '.  We  feel  that  any  successful  discriminant  should  not  depend  strongly  on 
the  local  geologic  setting  and  mining  practice.  In  paper  I  we  examined  events  that  occurred  in  central 
Asia,  in  this  chapter  we  consider  Scandinavian  events.  We  have  automated  the  algorithm  to  the  point 
where  discrimination  can  be  carried  out  solely  by  the  computer.  This  type  of  problem  has  also  been 
investigated  by  Baumgardt  and  Ziegler  (1989).  Their  approach  also  relics  heavily  on  the  expected  time- 
independence  of  spectral  modulation  in  the  coda  produced  by  ripplc-fired  events.  In  both  the  present 
work  and  Baumgardt  and  Ziegler  (1989)  the  underlying  prcmi.se  of  this  automation  is  that  in  the  future, 
if  lower  thresholds  arc  realized,  and  thus  a  significantly  greater  data  set  mast  be  examined,  it  will  be 
beneficial  and  desirable  to  distance  the  human  element  from  the  discrimination  process. 

1 .2  The  data  set 

The  bulk  of  the  data  ased  by  this  study  were  collected  by  sensors  in  the  NORESS  small- 
aperture  array  -  located  in  south-castem  Norway  -  from  1985  to  1986  (see  Figure  1.1  and  Table  1.1).  The 
NORESS  array  is  composed  primn.ily  of  25  vertical  component  sensors  deployed  roughly  2  m  deep  in 
shallow  vaults  arranged  in  a  set  of  concentric  rings  (Mykkeltvcit  et  al ,  1983).  The  fourth  and  outermost 
ring  is  roughly  3  km  across.  The  signal,  collected  by  GS-13  seismometers  which  have  a  flat  respon,sc  to 
ground  velocity  between  1  and  10  Hz.,  is  digitized  at  40  s  '.  NORESS  is  actually  part  of  a  significantly 
larger  array,  known  as  NORSAR,  and  is  situated  within  element  06C  at  this  array,  a  site  known  to  be 
particularly  sensitive  to  signals  propagating  from  Semipalatinsk  (Richards,  1988).  The  seismometers  are 
deployed  in  competent  igneous  rocks  of  granitic,  rhyolitic  and  gabbroic  composition  (Mykkeltvcit,  1987) 
and  Precambrian  or  Paleozoic  age  (Bungum  et  al .  1985).  The  site  is  thus  relatively  immune  to  the 
near-surface  resonance  of  seismic  energy.  A  more  complete  description  of  the  array  can  be  found  in 
Mykkeltvcit  et  al  ( 1983). 
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Event 

Latitude 

°N 

030 

59.31 

039 

59.31 

501 

58.34 

094 

59.73 

099 

61.55 

111 

60.19 

112 

58.90 

158 

58.34 

196 

58.34 

,522 

62.74 

523 

62.90 

208 

62.90 

216 

66.45 

524 

62.40 

525 

62.61 

504 

58..34 

505 

67.10 

506 

58.34 

526 

58.34 

239 

62.76 

247 

61.67 

266 

61.66 

270 

58.34 

298 

59.31 

507 

59.31 

508 

59.31 

509 

58..34 

510 

59.31 

511 

58.34 

512 

59.31 

513 

59.31 

514 

59.31 

515 

59.31 

516 

59.31 

517 

58.34 

518 

59.31 

519 

59.31 

520 

59.31 

503 

58.34 

521 

58..34 

407 

61.97 

472 

61.46 

Longitude 

Origin  Time 

Ml 

Event  type 

°E 

y/d-h:m:s  (UTC) 

6.95 

8.5/302-10:22:52.8 

1.9 

Blasjo  ex 

6.95 

8.5/3 10- 14:.50:5 1.4 

2.4 

Blasjo  ex 

6.43 

8.5/313-14:42:45.0 

— 

Titania  ex 

5.71 

85/331-04:53:32.1 

3.0 

earthquake 

4.65 

85/334-19:05:13.4 

3.0 

earthquake 

5.25 

85/341-14:15:43.2 

2.2 

earthquake 

5.98 

85/341-14:39:09.9 

1.9 

earthquake 

6.43 

85/365-13:36:49.6 

2.1 

Titania  ex 

6.43 

86/031-14:17:35.7 

1.9 

Titania  ex 

4.50 

86/036-23:35:41.0 

2.6 

earthquake 

4.86 

86/037-06:19:52.0 

2.3 

earthquake 

4.86 

86/037-06:20:05.4 

1.9 

earthquake 

14.89 

86/038-21:03:21.1 

2.2 

earthquake 

5.28 

86/044-13:39:00.0 

2.5 

earthquake 

5.07 

86/044-19:03:48.0 

2.6 

earthquake 

6.43 

86/045-14:13:25.0 

2.7 

Titania  ex 

20.60 

86/045-16:44:08.0 

2.6 

explosion 

6.43 

86/045-17:54:11.0 

2.3 

Titania  ex 

6.43 

86/047-18:19:41.0 

7.0 

earthquake 

5.29 

86A)57-02:11:.58  5 

1.9 

earthquake 

2.58 

86/067-16:21:18.3 

1.9 

earthquake 

4.53 

86A)89-03:22:48.7 

1.6 

carthqu:  ^<e 

6.43 

86/094-13:12:43.9 

1.9 

Titania  ex 

6.95 

86/120-10:18:48.2 

2.2 

Blasjo  ex 

6.95 

86/147-18:36:14.0 

2.3 

Blasjo  ex 

6.95 

86/148-17:51:57.0 

2.4 

Bla.sjo  ex 

6.43 

86/157-13:14:28.0 

1.7 

Titania  ex 

6.95 

86/170-03:55:08.0 

2.5 

Blasjo  ex 

6.43 

86/174-13.12:54.0 

1.8 

Titania  ex 

6.95 

86/191-20:10:42.0 

2.3 

Blasjo  ex 

6.95 

86/197-17:49:28.0 

2.3 

Blasjo  ex 

6.95 

86/204-^:47:10.0 

2.2 

Blasjo  ex 

6.95 

86/210-13:13:41.0 

2.3 

Bla.sjo  ex 

6.95 

86/2ll-17:59:.39.0 

2.4 

Blasjo  ex 

6.43 

86/226-13: 14:.39.0 

1.9 

Titania  ex 

6.95 

86/226-14-,39:57.0 

2.4 

Blasjo  ex 

6.95 

86/245-12:53:51.0 

2.1 

Blasjo  ex 

6.95 

86/252-17:55:58.0 

2.4 

Blasjo  ex 

6.43 

86/274-14:15:10.0 

1.9 

Titania  ex 

6.43 

86/282-14:13:52.0 

2.0 

Titania  ex 

2.33 

86/283  19:56:29. 1 

2.1 

earthquake 

.3.29 

86/299-11:44:54.1 

2.4 

earthquake 

Tabic  1.1.  Event  locations,  origin  timc.s,  local  magnitudes  and  t>'pes. 


In  audition  to  the  NORESS  data,  we  shall  use  an  event  recorded  by  the  NRDC  high-frequency 
stations  deployed  in  Kazakhstan.  U.S.S.R.  in  1987  (Given  et  al.,  1990).  The  recording  we  have  chosen 
is  of  the  calibration  explosion,  Chemex  2,  and  was  made  by  the  surface  sensor  at  Bayanaul. 

The  events  recorded  by  the  NORESS  array  consist  of  earthquakes  and  quarry  blasts  which, 
with  the  exception  of  one  event,  occurred  within  a  range  of  700  km  from  the  array.  Only  regional  events 
are  considered  here  sitKe  the  analysis  depends  on  the  retention  of  high-frequency  energy  in  the  coda. 
Alt  events  fall  within  a  local  magnitude  range  of  1.6  <  Ml  <  3.0.  Event  magnitudes,  locations,  origin 
times  and  idcntificatiorts  were  obtained  from  BaumgardI  and  Ziegler  (1988)  and  Screno  et  al.  (1987). 
All  frequency-spectral  estimates  have  been  computed  using  a  multi-taper  algorithm.  The  rationale  behind 
the  choice  of  this  algorithm  is  described  in  paper  1.  and  the  theory  describing  this  approach  can  be  found 
in  numerous  papers,  including  Park  et  al.  (1987)  and  Thompson  (1982). 

1.3  The  effect  of  tipple-firing 

At  least  at  the  macroscopic  scale,  the  practice  of  ripplc-hring  appears  to  have  Uttle  systematic 
effect  on  the  seismic  waveforms.  It  is  well  known,  however,  that  ripple-fired  events  tend  to  give  rise  to 
seismic  coda  possessing  highly  colored  spectra,  that  is,  spectra  enriched  in  power  in  certain,  preferred, 
frequency  bands  and  depleted  in  power  in  others  (Bell.  1977;  Baumgardt  and  Ziegler,  1988;  Stump  and 
Riemcr.  1988;  Smith  1989).  This  spectral  color  is  due  to  the  interaction  of  the  time-offset  wavefields  pro¬ 
duced  by  each  sub-explosion.  Brielly,  the  regular  repeution  and  superposition  of  similar  seismic  motions 
in  the  time  domain  leads  to  regular  amplification  and  suppression  of  power  in  the  frequency  domain.  The 
manner  in  which  the  wavefields  interact  undoubtedly  involves  nonlinear  processes;  however,  we  feel  that 
simple  linear  theory  is  sufficient  to  describe  the  most  obvious  result,  specifically  the  pronounced  spectral 
modulation.  As  dc.scribcd  in  paper  I.  and  by  numerous  other  authors  (Baumgardt  and  Ziegler,  1988; 
Sliuiip  and  Rcinke,  1988;  Smith,  1989;  Slump  ef  al ,  1989)  this  model  makes  the  a-ssumptions  that  the 
wavefields  produced  by  each  sub-explosion  w(t)  are  identical  and  that  they  suiKrpose  linearly.  Forcing 
all  shots  to  occur  at  regular  timc-micrvals  T  we  can  construct  the  wavelet  produced  by  the  ensemble  of 
sub-explosions  (lasting  a  total  of  D  seconds)  by  the  equation: 

r(t)  =  u(/)*[:i///(:5;)«(-^3)]  (1.1) 

where  *  represents  convolution.  Here,  III  is  the  shah  funetkm  (Bracewcll,  1986)  and  B  is  the  boxcar 
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function.  Hereafter  we  refer  to  this  representation  as  model  1.  By  Fourier  transforming  this  expression 
we  see  that  the  spectrum  of  the  entire  seismogram  equals  that  of  an  individual  event  multiplied  by  a  set 
of  cquispaced  sine  functions  -  collectively  referred  to  as  the  mtxlulation  function; 

X{f)  =  VV(/)[///(/7')  ,  (1.2) 

In  Figure  1 .2  we  display  the  modulation  function  resulting  when  39  sub-explosions  spaced  at  2S  ms  are 
superposed  in  this  manner.  Primary  reinforcement  occurs  at  multiples  of  40  H/.  (the  loci  of  the  main-lobes 
of  the  sine  functions).  The  side-lobes  have  insignificant  amplitudes  relative  to  the  main-lobes.  Tliey  can. 
however,  in  theory  allow  us  to  compute  the  duration  of  the  entire  quarry  blast.  The  duration,  D,  is  given 
by  the  inverse  of  the  width  of  a  single  sidclobc.  In  the  event  displayed  in  Figure  1.2,  this  value  is  0.975 
seconds,  the  known  duration  of  the  set  of  explosions.  As  discussed  in  paper  1.  using  the  model  described 
above,  we  predict  that  the  modulation  produced  by  ripple  firing  should  be  independent  of  time  in  the 
coda.  In  paper  1  wc  found  that  this  predicted  character  can  be  investigated  efficiently  by  the  computation 
of  frequency-time  displays  known  as  sonograms  (Markel  and  Gray.  1976;  paper  1).  In  Figures  1.3  and 
I A  arc  displayed  the  sonograms  computed  from  the  coda  generated  by  an  earthquake  and  a  quarry  blast 
respectively.  The  quarry  blast  (Figure  1.4)  clearly  shows  a  time-independent  spectral  modulation  whereas 
the  earthquake  (Figure  1.3)  does  not.  Often  the  two  types  of  events  do  not  contrast  as  well  as  these 
examples  do  when  presented  in  this  format.  For  this  reason  wc  have  found  it  bcnehcial  to  convert  the 
spectral  estimates  to  binary  form.  The  means  by  which  wc  accomplish  this  conversion  is  discussed 
fully  in  the  paper  1.  and  involves  comparing  a  relatively  unsmoothed  version  of  each  spectrum  with  a 
more  heavily  smoothed  one  that  resolves  only  the  large  scale  structure,  in  order  to  extract  the  regular 
modulation.  In  practice,  when  analyzing  the  events  considered  in  this  chapter,  we  simply  convolved 
the  spectra  with  boxcar  functions  .spanning  1 .0  and  2.5  Hz  respectively.  Wc  then  represent  all  sections 
of  the  spectra  where  the  UkuI  power  is  high  relative  to  the  more  regional  average  power  by  a  value 
of  -I- 1  and  where  it  is  low  by  a  value  of  -I.  In  this  manner  the  bulk  of  the  magnitude  information  is 
di.scardcd  and  the  spectra  arc  "llattcncd"  to  very  simple  binary  patterns.  When  analyzing  array  data, 
we  generalize  the  procedure  by  computing  such  a  binary  pattern  for  each  trace  individually,  and  then 
stacking  all  the  patterns.  Bccau.se  the  procedure  is  quite  nonlinear,  this  Ls  very  different  from  computing 
binary  sonograms  from  beams  as  in  Figures  1.3  and  1.4.  As  illusuatcd  below,  stacking  erfter  reduction 
to  binary  patterns  is  a  more  effective  approach  for  our  present  purpo.ses.  In  Figures  1.5  and  1.6  we 
display  array  stacks  of  the  binary  sonograms  computed  for  the  events  displayed  in  Figures  1 .3  and  1 .4 
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Figure  1.2  Spectral  modulation  predicted  for  an  event  consisting  of  39  sub-explosions 
located  at  the  same  point  in  space  and  offset  evenly  in  time  at  25  ms. 

respectively.  Since  typically  25  vertical  sensors  simultaneously  record  each  event,  the  values  in  these 
binary  stacks  typically  range  from  -25  to  +25.  The  original  spccu-al  estimates  have  been  corrected  for 
noi.se  by  subtracting  an  average  pre-event  sample.  Time -independent  spectral  modulation  is  present  after 
the  onset  in  the  coda  of  the  quarry  blast  only.  This  spectral  character  is  not  unique  to  this  event  but  is 
shared  by  virtually  all  the  events  identified  in  Table  1.1  as  explosions. 

1.4  The  cause  of  the  observed  spectral  modulation 

The  simplest  cxphination  of  the  observed  spectral  mixlulation  is,  as  discassed  in  the  previous 
section,  that  it  is  due  to  ripple-tiring.  The  main  lugumcnl  against  this  explanation  is  that  the  inferred 
delay  times  at  the  source  are  extremely  long.  Spectra  computed  from  a  typical  event  (030  -  pictured  in 
Figures  1  4  and  1 .6)  have  power  highs  spaced  at  roughly  5  H/  leading  to  an  inferred  average  shot  spacing 


0  10  20  30  40  50 

frequency  (Hz) 
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Figure  1.3  Seismogram  resulting  from  an  earthquake  kK'ated  342  km  from  NORESS 
(event  094)  and  corresponding  sonogram.  In  Figure.s  1.3  through  1.6  the  sonograms 
have  been  computed  from  a  stack  of  25  spectra,  each  computwl  from  an  individual 
vertical  channel  in  the  NORESS  array.  The  stacks  were  computed  after  offsetting  the 
seismograms  to  bcamform  for  the  Pp  phase.  In  addition  all  spectral  estimates  have  been 
corrected  for  noise  and  die  instrument  response.  The  spectral  amplitudes  in  figures  1.3 
and  1.4  arc  shown  on  a  linear  scale. 


of  200  ms.  In  paper  1  we  inferred  delay  times  as  high  as  400  ms  at  quarries  in  Karakhstan,  U.S.S.R. 
As  Baumgardt  and  Ziegler  (1988)  mention,  “slow  delays"  (from  500  to  1000  ms)  are  used  in  subsurface 
mining  where  the  intent  is  to  use  a  shot  to  remove  material  prior  to  the  next  shot.  We  have  reason  to 
believe,  however,  that  the  explosions  considered  in  paper  1  and  the  current  data  set  did  not  occur  in  the 
subsurface. 

With  the  aid  of  satellite  (SPOT)  photos  we  know  that  a  number  of  the  mines  in  Kazakhstan 
arc  at  the  free  surface  (Thurber  et  al,  1989).  The  Blasjo  explosions  arc  known  to  be  associated  with 
the  construction  of  a  d;im  (Baumgardt  :ind  Ziegler.  1988).  As  discus.scd  by  several  authors  (including 
luingefors  and  KihlstrOm.  1978)  the  short  delays  employed  at  frec-surfacc  mining  operations  generally 


9 


Figure  1.1  Seismogram  resulting  from  a  quarry  blast  l<x:alcd  301  km  from  NORES.S 
(event  030)  and  corresponding  sonogram. 

fall  in  the  range  from  1  to  100  ms  and  are  typically  on  the  order  of  20  to  30  ms.  Using  the  model 
descritxid  in  the  previous  section,  and  ,30  ms  off.scLs,  see  predict  spectral  amplification  at  multiples  of 
33  U/  -  well  beyond  the  Nyquist  frequency  of  the  NORHSS  data  set.  It  is  conceivable  that  the  closely 
spaced  mixlulations  (shown  in  Figures  1.4  ;uk1  1.6)  could  be  an  artifact  of  multiple-row  blasting  where 
short  delays  arc  used  between  successive  shots  in  each  row,  but  adj.iceMl  rows  arc  spaced  by  signilieanlly 
greater  delays.  Synthetic  experiments,  in  which  modulation  functions  arc  computed  for  a  variety  of 
quarry  blast  configurations,  suggest  that  this  is  a  plausible  argument;  however,  realistic  examples  taken 
from  the  literature  do  not.  For  example.  Stump  et  at  (1989)  describe  multiple-row  quarries  which  have 
interrow  time  spacings  of  42  ms.  This  argument  doc,s  not  rule  out  slow  delays,  either  between  successive 
shots  or  adj.icenl  rows,  but  suggests  that  we  should  look  for  altemalive  explanations  for  the  observed 
s()cctral  nuxlulalion. 

As  discussed  in  paper  1  and  Hedlin  et  at  (1488)  it  is  possible  for  a  wavciicid  to  acquire  a 


10 


lime  ajier  eveni 

Si)  iiM>  j:o  140  y<v7  iso 


Figure  1.5  S<  i>niojjram  a'suliiiig  Iront  ih>‘  (xirthgittke  prescnied  in  Figure  (0941 
and  ei>ires(KMiui(ig  binary  sonogram,  llte  ciMiversion  to  binary  foim  was  perfonncd  on 
each  channel  K'loir  slacking 

linie-independcni  s|vitral  inodulalii'n  during  (wopagalion  by  re.sonaiing  in  low  velocity  layers.  The  most 
likely  Kvaiions  oi  Liyer  reson.mce  ;ire  in  low  velocity  sediments  or  weathered  strata  near  the  ftee-surface 
close  to  the  source  and/or  the  array.  Considering  that  many  of  the  recorded  events  have  given  rise  to 
unmodulated  sjKClra  it  is  clear  that  no  signiticani  near  rcxciver  rc,sonance  is  liiking  place.  Furthermore, 
since  different  mcxiul.iiion  patterns  are  conmuMiiy  produced  by  different  events  with  the  same  location 
(such  as  successive  mine  explosions  at  the  s;une  mmol,  the  miHlulaiums  are  clearly  not  due  to  near  source 
resi'nance  We  cmicluik-  that  the  '[Vctral  m'ldulation  is  most  likely  due  to  intrinsic  source  processes. 

tinrd  explaiuition  relies  again  (  ii  source  multiplicity.  The  mixlulation  function  produced  by 
miHlel  I  IS  d‘>mu',i'.  (t  bv  tlK’  ir.  iin  lolx-s  oi  the  sm,  (unctions.  Tliese  are  the  only  features  tliat  can 
rc-alisiicalK  N-  ''[cs  i  d  lo  prixliicc  observ.ible  s(K'cir,il  peaks  when  the  time  delays  are  perfectly  reguhu'. 
MihIcI  I.  tiocci  vt  [  diH'  not  deserilx-  .1  verc  likely  c|u.inv  bl.isi.  As  discussed  by  many  authors  (including 
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Figure  1.6  Seismogram  resulting  from  the  quarry  blast  presented  in  Figure  1,4  (0.10) 
and  corresponding  binary  sonogram. 

Langefors  and  Kihlslrom,  1978:  Dowding,  198.*'  and  Stump  ct  ai.  1989)  ripple-fired  shots  in  quarries  are 
spatially  offset,  usually  in  a  regular  pattern.  At  each  shot  location  there  are  sometimes  several  vertically 
offset  (decked)  charge.s.  The  time-delays  between  the  .shots,  especially  in  multiple  row  blasting,  are 
not  necessarily  going  to  be  consistent.  Actual  shot  limes  often  deviate  a  considerable  amount  from  the 
intended  limes  (Slump  and  Reamer.  1988).  Kiuiwing  ific  near-surface  velixity  and  the  slowness  (p)  of 
the  energy  under  consideration  we  can  replace  actual  time  and  space  offsets  (iS/’,  and  ()A',,(^1',)  with 
apparent  time  offsets  by  noting  that. 

~  -  f')\rns0)  (*'/',  (i  ,i, 

The  luimuth  from  the  quarry  to  the  receiver  is  given  by  0.  All  ibe  aforementioned  factors  can  cause  a 
considerable  deviation  of  the  apparent  times  of  the  sub  explosions  fnim  a  common  value.  Using  these 
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apparcnl  time-offsets,  and  assuming  linear  superposition  and  commonality  of  basis  wavelet  w(t),  we  can 
construct  tlx;  wavelet  due  to  a  more  gencnil  quarrv  blast. 

n 

-  (1.4) 

I  -  I 

In  the  frequency  domain  this  expression  is  ajuivalcnl  to; 

n  »i 

^(/)  =  *  (53.so.(2rr/^7’“))'-*]^  (1.5) 

1-1  1-1 

Scatter  in  the  apparent  times  of  the  sub-explosions  reduces  the  dominance  of  the  main  lobes,  or  equiv¬ 
alently,  IcLs  the  side-lobes  ri.se  into  prominence  (paper  1).  To  illustrate  this  point  we  have  computed  a 
theoretical  modulation  pattern  for  a  quarry  blast  layout  adapted  from  that  of  a  real  life  quarry,  the  San 
Vel  quarry,  described  and  studied  by  Slump  and  Rcimcr  (1988)  and  Stump  et  al.  (1989).  As  displayed 
in  Figure  1 .7  the  sub-explosions  are  arranged  in  an  en  echelon  pattern.  The  shots  in  each  row  are  spaced 
at  25  ms  proceeding  from  west  to  east.  The  row  detonations  arc  separated  by  42  ms  proceeding  from 
south  to  north.  The  modulation  functions,  computed  for  energy  traveling  lo  observation  points  due  north 
and  east  of  the  quarry  with  a  slowness  of  1/7  s/km.  arc  displayed  in  Figure  1.8.  Although  the  dominant 
delay  time  is  25  ms,  Ihc  40  Hz  peak  docs  not  dominate  cither  modulation  function.  The  function  for  the 
station  to  the  north  can  be  constructed  by  multiplying  the  nuxlulation  function  due  to  13  shots  spaced  at 
25  ms  (representing  the  inicrshol  delays)  with  the  function  corresponding  to  3  shots  spaced  al  roughly  42 
ms  (representing  the  inter-row  delays  after  taking  into  account  the  delay  associated  with  the  propagation 
of  the  energy  between  the  rows).  The  two  functions  arc  in  competition,  and  the  result  is  that  the  broad 
main-lobes  of  the  latter  accentuate  the  side-lobes  of  the  former  lo  a  point  where  they  can  be  expected  to 
have  a  significant  impact  on  the  spectrum  of  the  quarry  blast. 

Using  a  technique  employed  in  paper  I  we  synthesize  a  quarry  blast  using  the  apparent  sub- 
explosion  times  occurring  in  the  event  described  above.  We  assume  a  common  waveform  is  generated 
by  each  sub-explosion  and  for  that  waveform  we  select  Ihc  calibration  explosion  Chemex  2  detonated  in 
Kazakhstan.  U.S.S.R.  and  recorded  at  the  station  at  Bayanaui.  (We  have  resorted  to  this  data  set  simply 
because  the  40  llz  NORESS  data  do  not  have  adequate  resolution  in  time  lo  permit  the  millisecond 
offsets  required  by  this  quarry.)  The  Chemex  2  recording  was  made  al  250  s  '.  This  “Green’s  function” 
is  linearly  stacked  upon  itself  3*^  limes  after  including  the  offsets  appropriate  for  the  observation  point 
due  north  of  the  quarry.  Although  we  have  chosen  to  create  Ihc  synthetic  quarry  blast  by  offsetting  and 
stacking  a  Green's  function  in  the  time  domain,  the  equivalent  result  could  be  achieved  by  multiplying 
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Figaro  1.7  The  layout  of  sub-explosions  in  an  en  echelon  quarry  blast.  Shooting 
within  the  rows  is  spaced  in  time  at  25  ms.  Adjacent  rows  arc  separated  by  42  ms  in 
lime.  This  pattern  is  adapted  from  Stump  ei  al.  (1989). 


the  .spectrum  of  the  Green’s  function  by  the  complex  modulation  function  which  underlies  the  solid  curve 
pictured  in  Figure  1.8.  Prior  to  computing  the  .sonogram,  tlie  "synthetic"  seismogram  was  low-pass  filtered 
between  0  and  20  Hz  and  decimated  to  one  point  in  5  to  mimic  a  NORESS  recording.  The  sonogram 
(displayed  in  Figure  1.9)  is  dominated  by  lime- independent  structure.  Assuming  this  modulation  pattern 
was  due  to  main-lobe  activity,  one  would  estimate  a  dominant  delay  time  to  be  roughly  170  ms  (the 
inverse  of  6  Hz).  We  know,  however,  that  this  structure  is  due  to  side-lobe  activity  and  is  controlled 
in  thi.s  case  by  the  total  duration  of  the  quarry  blast.  Because  of  the  manner  in  which  the  sonogram  is 
calculated,  the  frequency  estimates  arc  heavily  smixithcd.  Longer  time  windows  would  allow  a  more 
accurate  estimate  of  the  frequency  spiicing  of  the  modulation.  In  fact  we  know  that  the  average  sidclobc 
widlh  is  roughly  2.6  Hz  (sec  Figure  1.8)  and  (hat  the  duration  of  the  quarry  blast  is  384  ms. 

It  seems  that  there  is  a  fundamental  ambiguity  in  the  spectral  modulation  produced  by  ripple- 
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Figure  1.8  The  amplitude  of  modulation  functions  resulting  from  the  shot  introduced 
in  Figure  7.  The  solid  and  dashed  curves  represent  energy  traveling  at  a  slowness  of 
1/7  s  km  'to  stations  due  north  and  east  of  the  quarry  respectively. 


fired,  and  hence  non-instantaneous.  events.  Without  a  priori  information  about  what  occurred  at  the 
source  we  cannot  be  sure  if  the  modulation  spacing  is  controlled  by  the  duration  of  the  entire  set  or  by 
the  dominant  inter-shot  apparent  time  spacing.  This  experiment  shows  that  our  discriminant,  perhaps, 
will  not  recognize  quarry  blasts  because  they  arc  ripple  fired  per  sc.  but  because  they  last  an  intermediate 
length  of  time.  Instantaneous  events  give  rise  to  unmodulated  spectra.  Extremely  long  events  (for 
example  large  earthquakes)  should  produce  very  finely  modulated  spectra,  such  that  the  modulation  is 
masked  by  scattering  and  noise. 

In  paper  1  we  considered  two  types  of  events  -  calibration  explosions  which  were  detonated 
by  American  and  Soviet  scientists  (Given  et  at ,  1990)  and  did  not  involve  ripple  firing.  Using  a  priori 
information  we  strongly  suspected  the  rest  of  the  events  were  quarry  blasts.  This  information  included 
satellite  photos,  provided  by  Prof.  Clifford  Thurber  (at  the  University  of  Wisconsin),  which  showed 
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Figuro  1.9  Time  scries  and  sonogram  representing  a  "synthetic"  quarry  blast.  The 
synthetic  was  constructed  by  linearly  stacking  a  sei.smogram  produced  by  the  Chemex 
2  explosion  in  KazakhsUtn  U.S.S.R.  upon  itself  after  offsetting  in  lime  to  mimic  the 
quarry  pictured  in  Figure  1.7  w'hcn  observed  from  a  point  to  the  north  (see  Figure  1.8). 
TTie  original  seismogram  was  sampled  at  250  s'.  The  synthetic  was  low  pass  filtered 
and  decimated  to  simulate  recording  conditions  similar  to  the  NORESS  array. 


surhite  mining  activity  in  the  vicinity  of  some  of  the  events.  In  addition  it  is  known  that  the  region  has 
iillle  natural  seismicity  (Leith.  1987).  Time  independent  special  modulation  was  only  observed  in  the 
latter  set  of  events  and  was  attributed  to  the  source  multiplicity.  The  present  study  and  the  previous  one 
are  consistent  in  suggesting  that  quarry  blasts  can  be  discriminated  from  non  rippIc-fircd  events. 


1..5  The  automatic  discriminant 


For  our  purposes  it  is  irrclev;uil  whether  the  time  independent  spectral  features  observed  in  the 
coda  produced  by  quarry  blasts  are  due  to  main  lobes  or  to  sidclobcs  in  the  modulation  functions.  Ripple 
tired  events  tend  to  give  rise  to  lime  independent  spectral  modulation,  the  earthquakes  examined  in  this 


study  do  not.  To  examine  this  modulation  we  have  developed  a  means  to  “expand”  a  time  series  into 
a  matrix  of  numbers  depending  on  frequency  and  time.  Typical  patterns  obtained  from  recordings  of 
an  earthquake  and  a  quarry  blast  (Figures  1.5  and  1.6)  illustrate  that  it  can  be  very  easy  to  discriminate 
visually  between  these  two  types  of  events  given  these  time-frequency  displays.  In  paper  1,  asing  the  same 
approach,  we  found  a  similar  degree  of  success  in  discriminating  between  quarry  blasts  aiKl  single-event 
explosions.  Given  the  current  interest  in  the  problem  of  discriminating  quarry  blasts  from  earthquakes 
and  single-event  explosions  and  the  large  numbers  of  events  involved,  we  feel  it  is  important  to  extend 
the  algorithm  so  that  human  intervention  is  distanced  from  the  discrimination  process,  to  a  point  where 
the  patterns  can  be  rccognb.cd  automatically  by  the  computer.  One  method  we  have  found  to  be  very 
effective  involves  the  computation  of  a  two-dimensional  Fourier  transform  of  the  sonogram  matrices. 
This  can  be  considered  as  an  extension  of  cepstral  analysis  (Tribolet,  1979).  In  the  standard  cepstral 
analysis  a  Fourier  transform  of  the  log  of  the  amplitude  spectrum  is  computed  to  highlight  any  regular 
spectral  modulation  regardless  of  its  longevity.  The  independent  variable  is  Imown  as  the  quefrency  and 
has  units  of  time.  The  form  of  cepstral  analysis  we  are  proposing  is  more  demanding,  however.  A 
given  point  in  the  2-D  cepstral  matrix  not  only  represents  spectral  modulation  at  a  certain  quefrency.  but 
periodic  along  the  time  axis  at  a  certain  frequency.  It  is  thus  a  simple  matter  to  isolate  energy  periodic 
in  frequency  and  independent  of  time. 

To  illastrate  our  point  we  display  two  2-D  cepstra  in  Figures  1.10  and  1.11.  The  first  was 
computed  from  the  first  100  seconds  of  coda  of  event  030  (Figures  1.4  and  1.6).  the  second  was  computed 
from  the  coda  of  the  earthquake  094  (Figures  1.3  and  1.5).  The  quarry  blast  has  significantly  more 
energy  at  zero  frequency  (along  the  time  axis)  than  the  earthquake.  The  quefrency  at  which  the  power  is 
concentrated  in  the  2-D  cepstrum  computed  from  the  coda  produced  by  the  quarry  blast  is  roughly  0.2 
.seconds  (reflecting  the  spectral  modulation  with  5  Hz  spacing.  Slices  at  zero  lime-frequency  through  2-D 
cepstra  computed  from  the  coda  produced  by  a  quarry  blast  (event  .507)  and  all  the  earthquakes  in  the 
data  .set  are  shown  in  Figure  1.12.  As  expected,  the  quarry  blast  is  a  singular  event.  The  most  nobceable 
feature  in  the  quarry  blast  cepstrum  is.  obviously,  the  significant  peak  at  a  quefrency  of  0.2  seconds. 
We  expect  that  ripple-fired  events  should  give  ri.se  to  significantly  larger  extreme  cepstral  values  than 
earthquakes.  Supporting  this  thesis  are  the  histograms  in  Figure  1.13  showing  the  observed  distributions 
of  cepstral  extremes  for  the  entire  earthquake  and  quarry  blast  populations  examined  in  this  study. 

Although  we  are  most  interested  in  the  quarry  blast  cepstra.  wc  can  gain  some  important 
insight  from  the  earthquake  cepstra  which  illustrate  the  2-D  cepstral  structure  that  can  be  expected  in  the 
absence  of  .source  multiplicity.  These  cepstra  show  what  time-independent  structure  will  be  acquired  by 
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Figure  1.10  Two-dime  as  ion  al  cepstnim  computed  from  the  coda  produced  by  a  quarry 
blast  (event  0,10).  The  first  100  seconds  of  the  coda  were  considered. 

a  propagaiinj;  wavelet  or.  in  otlter  words,  they  arc  indicative  of  the  region's  natural  level  of  resonance. 
We  propirse  to  identify  events  as  quarry  blasts  by  searching  for  anomaloasly  high  global  extrema  in  the 
time  independent  segments  of  the  2  D  cepsira.  To  calibrate  the  algonthm.  to  account  for  the  natural 
resonance  in  the  region,  we  make  the  judgment  of  what  is  a  large  value  on  the  basis  of  what  extrema 
non-rtpple-fircd  events  produce.  The  consideration  of  global  extrema  in  these  2-D  cepstra  is  a  problem 
that  is  well  suited  for  analysis  using  the  statistics  of  extremes  (Gumbcl.  19‘i8). 

In  Figure  1.1.1  it  is  clear  that  the  logs  of  the  extreme  amplitudes  are  centrally  distributed,  and 
there  are  no  significant  outliers.  The  Kolmogorov-Smimov  test  suggests  that  tlic  earthquake  cc-pstral 
extremes  follow  a  log-normal  distribution.  However  we  would  like  to  avoid  the  adoption  of  a  specific 
underlying  distribution  since  we  have  no  fundamental  reason  forchcxising  one  and  since  we  only  have  16 
earthquakes.  It  is  known  (e.g..  Gumbcl.  19.S8  and  Wei.ssman.  1978)  that  when  dealing  with  observations 
of  extreme  values  the  underlying  distribution  need  not  be  assumed,  but  the  behavior  can  be  modeled  using 
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Figure  1,11  Two-dimensional  cepstnim  computed  from  the  coda  produced  by  an 
earthquake  (ev’nt  094).  The  first  100  .seconds  of  the  ctxJa  were  considered. 

functions  that  are  asymptotically  valid  as  the  number  of  .samples  examined  and  the  number  of  points  in 
each  sample  approaches  infinity  (Kennedy  and  Neville.  1974).  Selecting  the  exponential  asymptote,  the 
cumulative  probability  (P)  that  an  extremum  belonging  to  the  earthquake  population  will  be  less  than  the 
one  observed  is  given  1. . 


P  =  (1.6) 

where  the  expression  for  the  reduced  viiriate  (y)  is: 

y  —  a(u  —  u)  (1-7) 

The  terms  a  and  u  are  the  dispersion  parameter  and  the  mode  of  the  distribution  respectively  and  are 
estimated  directly  from  the  population  of  earthquake  extremes  sltown  in  Figure  1.13  (Kcnncrly  and 
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Figiirr  1.12  Slices  through  (he  time-independent  portions  of  two-dimensional  cepstra 
computed  from  a  quarry  blast  (event  S07)  and  all  the  earthquakes  considered  in  this 
study.  The  quarry  blast  is  shown  as  the  solid  line. 


Neville.  1974).  The  log  of  the  ccpstral  extremum  of  interest  is  represented  by  u.  We  Qnd  that  for  this 
distnhulion  a  and  u  equal  7.204  and  3.55,  respectively. 

Given  this  probability  function  we  can  pose  the  discrimination  problem  in  terms  of  a  standard 
hypothesis  test:  Let  the  null  hypothesis  (//o)  be  that  a  newly  recorded  event  belongs  to  the  population  of 
earthquakes  used  to  calibrate  the  technique.  If  the  ccpstral  cxU’cmum  calculated  for  this  event  exceeds 
a  certain  threshold  determined  fnim  the  distribution  (1.6),  then  we  can  reject  the  null  hypothesis  (Hq) 
at  a  preset  confidence  level,  and  conclude  that  the  event  is  probably  a  quarry  blast.  For  example,  on 
Figure  1.14.  this  threshold  was  selected  such  that  for  points  that  plot  above  the  threshold  line,  the  null 
hypothesis  is  rejected  with  only  a  5‘^  risk  of  doing  so  erroneously.  In  other  words,  we  state  that  events 
above  the  line  do  not  belong  to  the  earthquake  population,  al  the  95%  confidence  level.  In  spite  of  the 
apparent  c'^ticiency  of  the  discriminant  illustrated  on  Figure  1.14,  we  must  remember  that  the  calibration 
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Figure  1.13  Histograms  showing  the  observed  distributions  of  global  extreme  cep- 
stral  values  computed  from  the  cc^  produced  by  earthquakes  (top)  and  quarry  blasts 
(bottom). 

of  the  distribution  (8)  is  based  on  our  (small)  sample  of  16  identified  earthquakes,  so  that  the  test  is  in 
fact  “data  fitted”.  Confirmation  of  our  claim  of  success  will  have  to  he  based  on  an  independent  sample. 
In  this  figure  the  symbol  size  is  directly  proportional  to  the  signal  to  noise  ratio  (derived  from  average 
spectra  encompassing  the  time  from  50  seconds  before  and  after  the  compressional  onset).  Of  the  26 
quarry  blasts  considered,  23  lie  above  the  95%  confidence  level.  Of  the  two  that  fall  well  below  this 
limit,  one  (event  .505  -  located  in  northern  Sweden)  had  extremely  low  signal  to  noi.se  levels  (less  than 
10  dB)  and  the  other  (event  504)  produced  only  a  very  broad  specual  modulation.  The  three  earthquakes 
located  above  a  probability  of  0.8  (events  1 1 2,  523  and  208)  all  suffered  from  signal  to  noise  ratios  less 
than  10  dB. 
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Figure  1.14  The  cumulative  probabilities  of  extreme  cepstral  values  derived  from  the 
coda  produced  by  all  events  in  the  data  set.  The  quapy  blasts  are  denoted  by  octagons, 
the  earthquakes  are  represented  by  stars.  The  likelihood  that  the  assumption  that  the 
event  is  an  earthquake  is  invalid  increases  with  this  probability.  For  points  above  the 
0.95  threshold,  the  hypothesis  that  the  corresponding  events  belong  to  the  earthquake 
population  is  rejected  at  the  5%  risk  level.  The  event  number  (along  Ihe  horizontal 
axis)  indicates  Ihe  location  of  the  event  in  table  1.1.  The  symbol  size  scales  with  the 
signal  to  noise  ratio  (see  insert). 


1.6  Conclusions 

In  paper  1  we  advanced  the  preliminary  observation  that  ripple-fired  events  lend  to  give  rise 
to  coda  dominated  by  time- independent  spectral  features  and  that  thi.s  quality  should  be  exploited  to 
discriminate  these  events  from  earthquakes  and  single -event  explosions. 

In  this  chapter  we  have  demonstrated  that  this  can  also  be  done  with  a  high  degree  of  success 
when  considering  earthquakes  .ind  quarry  blasts.  We  have  found  that  quarry  blasts  tend  to  produce 
modulated  spectra,  but  the  modulations  may  not  result  directly  from  the  ripple-firing;  they  may  exist 
simply  because  the  event  is  non-instantaneous.  We  have  produced  an  empirical,  calibrated,  approach  to 
the  discrimination  problem  which  allows  for  local  seismic  resonance.  We  have  automated  the  approach 
to  tlic  point  where  discrimination  can  be  carried  out  solely  by  the  computer.  We  have  examined  a  data 
set  consisting  of  26  quarry  blasts  and  16  earthquakes  and  have  found  that  with  few  exceptions  the  two 
populations  are  well  separated  by  our  approach.  The  events  which  failed  to  be  identified  with  a  high 
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degree  of  confidence  generally  suffered  from  low  signal  to  noise  ratios. 

By  comparing  our  current  results  with  those  in  the  earlier  work  we  have  illustrated  the  ability  of 
the  algorithm  to  accommodate  changes  in  the  recording  env'u^onmcnt.  local  geologic  setting  and  mining 
practice.  Based  on  the  results  presented  in  paper  1,  we  expect  that  we  would  have  a  similar  degree  of 
success  in  discriminating  between  single-event  explosions  and  quarry  blasts. 
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Chapter  2 


The  migration  of  small-aperture  array  recordings 
to  image  local  crustal  scatterers 

Abstract 

This  chapter  is  concerned  with  the  development  of  a  technique  to  scan  seismic  coda  recorded  by 
single  small-apeiture  seismic  arrays  for  phases  generated  locally  by  scattering  from  large  heterogeneities, 
or  topographic  undulations.  In  essence,  a  hyperboloid-summation  lime  migration  is  performed  on  array 
raordings  of  distant  primary  events  to  enhance  local  scattered  pha,ses.  This  technique  can  be  applied  to 
single  events  but  should  give  more  robust  results  when  applied  to  a  broadly  distributed  suite  of  primary 
sources.  Prior  to  analyzing  recorded  data,  synthetics  arc  examined  to  assess  the  spatial  resolution  that  can 
be  achieved  with  the  imaging  algorithm.  A  preliminary  analysis  of  real  data  is  performed  to  provide  initial 
images  of  the  distribution  of  scatterers  in  the  vicinity  of  the  NORESS  array  in  southern  Norway.  Several 
telescismic  events  are  analyzed  separately  to  assess  the  feasibility  of  using  a  widely  distributed  suite  of 
dissimilar  primary  sources  in  the  imaging  process.  The  analysis  indicates  that,  due  to  the  protracted  nature 
of  most  telcseisms.  deconvolution  of  the  event  records  should  be  done  before  a  final,  spatial,  image  of  the 
scaitca-rs  can  be  produced.  The  initial  analysis  of  recorded  data  indicates  that  stable  apparent  secondary 
seismic  sources  are  present  in  the  vicinity  of  the  NORESS  array.  Thc.se  are  tentatively  interpreted  as 
scatterers  excited  by  the  primary  events. 


2.1  Inuoduction 

One  of  the  important  aspects  of  .seismic  monitoring  is  to  understand  the  generation  of  seismic 
coda,  particularly  the  near-station  mechanisms.  Two  mechanisms  that  have  been  investigated  include 
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seismic  resoiuncc  (in  low  vcUxrity  strata)  and  scattering  by  inhumogeneitics  and  topographic  undulations 
(both  at  the  free  surface  and  at  buried  interfaces).  The  two  pnxresses  appear  to  have  different  relative 
importance  in  the  continental  and  oceanic  crusts.  Recent  w(x^lc  (by  Screno  and  Orcutt,  1985a.b;  Screno 
and  Orcutt,  1987;  Mallick  and  Fra/er.  1990)  has  produced  evidence  that  in  the  oceanic  crust,  coda  waves 
arc  likely  dominated  by  rcsorrancc  in  the  water  and  sedimentary  horizons.  Other  studies  (eg.  Aki,  1969; 
Aki  and  Chouct,  197.^)  have  argued  that  coda  waves  generated  in  the  continental  crust  are  most  likely 
due  to  scattering  by  heterogeneities.  Although  studies  quantifying  the  effects  of  scattering  have  mostly 
been  statistical  and  have  dealt  with  the  influence  of  small-scale  random  scaltcrers.  some  deterministic 
studies  (Key,  1967;  Key,  1968;  Gupta  et  al .  I990a,b;  Bannister  et  al .  1990;  Lay,  1987;  Lynnes  and 
Lay,  1989)  have  produced  compelling  evidence  that  large  features  arc  capable  of  producing  significant 
amounts  of  scattered  ettergy,  principally  in  the  ftxm  of  large,  idcntihabic  seismic  phases,  which  affect 
nuclear  monitoring  and  discrimination:  These  prominent  scattered  phases  may  be  confused  with  direct 
arrivals.  Being  able  to  identify  these  'secondary"  sources  is  the  first  step  toward  being  able  to  reduce 
this  confusion. 

Several  attempts  have  been  made  to  l(x;ate  the  sources  of  scattered  phases  both  near  the  receiver 
(Key.  1967;  Gupta  ft  al .  1990a.h;  Bannister  et  al ,  1990)  and  near  the  source  (Lay,  1987;  Lynnes  and  Lay. 
1989).  The  work  to  dale  suggests  that  the  most  significant  sources  of  identifiable  phases  are  topographic 
features  at  the  free-surface. 


2.2  The  imaging  technique 

The  problem  of  imaging  sources  of  scattered  energy  using  two-dimensional  (2-D)  array  records 
bears  a  strong  rc.scmblancc  to  the  use  of  1-D  seismic  reflection  profiles  to  produce  images  of  subsurface 
velocity  contrasts.  As  a  result,  we  have  drawn  on  the  experience  of  the  reflection  seismology  community 
in  the  current  analysis. 

Recently,  we  (Hedlin  et  al ,  1990)  have  examined  a  suite  of  synthetic  and  recorded  events 
(observed  by  the  NORESS  array  in  Norway)  and  developed  a  systematic  technique  to  image  nearby 
large  scaltcrers.  The  underlying  premise  is  that  incident  seismic  waves  impinging  on  a  scattcrer  generate 
coda  waves  recorded  by  the  array.  From  this  point  of  view  each  scattcrer  may  be  treated  as  a  secondary 
source  which  is  excited  with  a  delay  estimated  from  elementary  ray  theory.  We  have  applied  a  modified 
bcamforming  technique  to  array  records  to  enhance  signals  radiated  by  such  faint  secondary  sources. 
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2  2  .0  The  adaptation  of  hyperbola  summation  migration  to  small  aperture  array  data 

To  describe  the  technique,  wc  adopt,  for  the  time  being,  a  simple  model  describing  the  origin 
of  locally  scattered  waves  -  one  in  which  a  single  omnidirectional  point  scattcrer  exists,  is  impulsively 
excited  at  a  time  t=0  and  produces  seismic  motions  that  arc  recorded  by  aii  array  of  sensors  located 
at  the  free  surface  (z=0;  x=Xj,  y=yj)  where  j  varies  from  1  to  N  (the  number  of  sensors  in  the  array). 
Furthermore,  we  assume  that  the  scattered  energy  propagates  at  a  constant  velocity.  V,  and  there  is 
no  dispersion.  Without  loss  of  generality  we  place  the  scattcrer  at  x=0,y=0  and  find,  using  simple 
trigonometry,  that  the  time  tj  at  which  the  scattered  wavefront  should  pass  the  jth  sensor  is  given  by: 


(2.1) 


or,  defining  -  ( ^ ).  the  vertical  travel  time  from  the  scattcrer  to  the  free  surface,  equation  2.1  becomes: 

(2.2) 
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Clearly,  the  surtace  in  ^  dimensional  space-time  which  describes  the  arrival  limes  predicted  for  energy 
arriving  from  the  .scattcrer  is  a  hyperboloid  of  revolution  (whose  axis  of  symmetry  is  the  lime  axis). 
From  examination  of  equation  2.1  it  is  made  clear  that  when  the  scattcrer  is  located  in  the  x-y  plane  the 
hyperboloid  degrades  to  a  circular  cone  with  its  apex  located  at  the  origin.  If  the  space-lime  location  of 
this  sensor  is  oubidc  the  hyperboloid,  the  sensor  cannot  have  been  influenced  by  the  scattcrer.  Conversely, 
if  the  sensor  is  located  within  the  hyperboloid,  the  scattered  wavcficld  will  have  already  passed.  Clearly, 
given  this  simple  situation,  to  achieve  the  greatest  enhancement  of  seismic  motions  caused  by  this  source 
(at  the  expense  of  motions  caused  by  secondary  sources  at  other  locations),  and  thus  to  achieve  the 
Kst  image  of  the  source  itself,  one  should  sum  the  motions  recorded  by  the  sensors  at  times  that  will 
reposition  each  .sensor  in  .space -time  onto  the  surface  of  the  hvpcrtxiloid.  This  is  a  simple  extension 
of  the  hyperbola  summation  migration  method  (Yilmaz,  1987)  in  seismic  reflection.  Since  a  delay-and- 
sum  operation  is  involved,  this  mcthixl  is  akin  to  standard  beamforming  techniques  and  bears  some 
resemblance  to  f-k  analysis.  This  technique  is  distinguished,  however,  by  its  cognizance  of  not  only  the 
wavenumber  but  the  onset  time  of  a  seismic  arrival.  As  a  result,  the  technique  can  be  u.sed  to  infer  the 
likely  geographic  liKation  of  the  seconcLu^  .source.  This  mcthixJ  differs  from  previous  atiempls  to  locale 
large  near-recciver  scatlerers  since  it  is  capable  of  acccwnnuxlating  simultaneously  many  primary  events 
from  Jiflerent  azimuths,  to  give  a  balanced,  redundant,  illumination  of  local  scatlerers  while  suppressing 
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the  influence  of  near-source  scattcrers. 

The  method  has  been  developed  to  accommodate  secondary  sources  located  at  any  point  in  the 
crust  but,  in  view  of  the  findings  made  by  previous  studies  (Gupta  et  al.,  1990a,b;  Bannister  et  al.  1990), 
the  following  analysis  is  directed  toward  sources  of  R,  energy  located  at  the  free  surface 


2.2.b  Processiriff  details 

To  image  a  region  of  the  crust,  we  subdivide  the  area  into  small  segments  and  considers 
them  individually,  in  sequence.  It  is  assumed  that  a  secondary  source  exists  within  the  sub-area  of 
current  interest.  Since  the  location  of  the  hypothetical  scatterer  is  known,  the  adoption  of  a  slowness  of 
propagation  of  the  scattered  waves,  p^,  allows  the  space-time  hyperboloid  to  be  defined.  For  each  crustal 
sub-area  being  scanned  ‘hf  ecordings  of  the  teleseismic  event  are  summed  after  shifting  them  in  time  as 
required  by  the  hvi'c'  <  .oid.  In  practice,  since  it  is  not  known  when,  in  absolute  time,  excitation  occurs, 
the  time  axis  is  translated  .so  that  the  onset  time  of  the  primary  energy  at  the  array  occurs  at  t=0.  Using 
clemenL'sy  ray  theory  it  is  possible  to  estimate  what  time  delay,  t,  should  exist  between  the  arrival  of  the 
prirnary  energy  and  the  .scattered  phases  originating  al  the  hypothetical  scatterer  after  excitation  by  the 
primary  energy.  As  is  illustrated  in  Figure  2.1,  considering  a  single  event-scatterer  pair,  the  time  offset, 
r,  between  the  arrivals  of  energy  propagating  directly  from  the  source  (at  vector  slowness  p,)  and  via 
the  scatterer  (vector  slowness  p,)  at  a  vector  distance  R„,  from  the  array  is  given  by; 


r  =  Ra.  "(p. -p.,)  (2.3) 

It  is  possible  to  estimate  p,  by  coasidering  a  suite  of  broadly  distributed  events  and  computing  a  number 
of  preliminary  images  while  slowly  varying  this  parameter  and  selecting  the  value  that  brings  the  image 
into  the  sharpest  focus.  The  .slowne,ss  of  the  incident  energy,  p,.  is  well  constrained  by  fitting  a  least- 
squares  best  fit  plane  (in  x-y-l  space)  to  the  first  breaks  on  the  array  records.  By  systematically  scanning 
the  crustal  volume  about  the  array,  we  can  generate  an  image  which  Is  interpreted  as  a  map  of  local 
scattcrers.  To  date,  we  have  considered  scattering  interactions  excited  by  P  waves,  using  both  synthetic 
seismograms  and  recorded  data. 

In  practice,  there  arc  a  number  of  complicating  factors  that  must  be  taken  into  account  prior  to 
constructing  images  of  faint,  local  sources  using  synthetics  or  recorded  data. 

I)  If  the  .scattered  field  coasisls  of  surface  waves,  propagating  al  the  surface  of  a  crust  whose 
velocity  increa.se.s  with  depth,  no  single  velocity,  V.  will  allow  a  comprehensive  description  of 
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Figure  2.1  The  vectors  appearing  in  equation  2.1  are  illustrated  in  this  figure.  The 
incident  wavefront  arrives  from  the  north  at  a  vector  slowness  of  p,.  A  scattered 
s^aselield  originates  within  the  small  box  (3  km  on  a  side)  and  propagates  outward  to 
the  NORESS  array  (dark  circle  at  the  center  of  the  figure)  at  vector  slowness  p,.  The 
vector  from  the  array  to  a  sample  block  is  Ro,.  Both  the  array  and  the  sample  block 
arc  shown  to  scale. 


the  scattered  wavefield  -  the  energy  will  become  dispersed.  As  a  result  it  becomes  necessary 
to  integrate  over  a  span  of  lime  centered  on  the  space-time  hycrboloid  which  is  defined  by  a 
velocity  representative  of  the  average  speed  of  propagation  of  the  surface  wave  energy.  We 
have  used  integration  time  windows  of  up  to  5  seconds  in  duration. 

2)  Based  on  the  results  of  synthetic  tests  we  have  not  attempted  to  use  the  full  coherent,  or 
phased,  recordings.  Because  we  arc  arc  dealing  with  dispersive  energy,  we  have  concluded 
th.it  the  trade-off  between  image  resolution  and  stability  can  best  be  mitigated  by  migrating 
coherent  records  hut  convening  the  results  to  envelopes  prior  to  integration  in  a  method  akin 
to  incoherent  bcamforming  (Ringdal  el  al .  19751. 

3)  When  considering  the  recorded  data,  it  is  immediately  apparent  that  local  scattered  phases. 


excited  by  tcleseismic  energy,  canned  be  examined  in  the  absence  of  this  (primary)  energy.  It 
is  necessary  to  mask  the  more  energetic  primary  sinircc  by  suppressing  the  energy  that  has 
propagated  directly  fnrm  it.  The  method  wc  have  adopted  to  accomplish  this,  intnxluced  by 
Gupta  el  (j/  (P>9()a),  is  known  as  beam  coirtx'lion  ;uid  simply  involves  the  coherent  subtntclion 
of  the  primary  source  beam  from  all  the  individual  channels  to  yield  residual  sei.smograms. 

4)  The  recorded  signals  are  immersed  in  seismic  noise,  much  of  which  is  microseismic  energy 
originating  in  the  North  Sea  (Fycn,  1986),  and  thus  exists  in  a  relatively  low  frequency  band 

to  Hz).  Fortunately,  the  surface  wave  energy  of  interest  exists,  mainly,  in  the  band  from 
to  3.0  Hz,  ;tnd  thus  the  signal  to  noi.se  ratio  (SNR)  can  he  significantly  improved  by  narrow 
band  pa.ss  filtering  of  the  array  records. 

5)  Since  wc  arc  considering  surface  wave  propagation,  a  gain  proportional  to  y^|R„|  should  be 
applied  to  compensate  for  amplitude  loss  due  to  cylindrical  geometrical  spreading. 

6)  Finally,  further  complication,  that  exists  only  in  the  recorded  data,  is  caused  by  the  non- 

t 

impulsive  nature  of  most  tclescisms  -  they  consist  of  a  protr;ictcd  sequence  of  arrivals.  Discus¬ 
sion  of  the  impact  of  this  complication  is  deferred  to  scclitMis  2.3  and  2.4. 

2.3  Analy.si.s  of  synthetic  seismograms  -  imaging  resolution 

The  resolution  of  the  imaging  technique  can  be  gauged  from  experiments  with  synthetics. 
Using  wavenumber  integration  (ApscI  and  Luco.  1983;  Luco  and  Apscl,  1983)  we  consider  a  localized 
omnidirectional  scattcrer  (delta  function  in  space)  illuminated  by  an  impulsive  incident  wave  (delta 
function  in  time),  and  compute  synthetic  vertical  seismograms  individually  for  each  of  the  sensors  in 
the  array.  Since  no  energy  from  the  primary  source  is  added  to  the  local  synthetic  phases,  this  yields  a 
perfectly  beam -corrected  image  (Figure  2.2).  This,  and  all  subsequent  images,  have  been  computed  by 
subdividing  the  area  shown  into  1681  bhKks  3  km  on  a  side.  In  the  preliminary  analysis  of  synthetic 
sei.smograms  ;ind  recorded  data  we  have  chosen  a  surface  wave  slowness  of  .304  s/km  (the  slowness  of 
the  dominant  phase  in  the  synthetic  Ry  packet)  .and  have  bamlpass  filtered  the  records  between  1  and  3 
Hz.  All  synthetic  images  have  been  adjusted  for  the  amplitude  loss  due  to  gceimctrical  spreading. 

As  illustrated  in  Figure  2.2  ttie  r;Klial  resolution  of  the  point  source  Ls  limited.  This  is  due  to 
some  extent  to  the  dispersion  cif  die  surfiice  wave  ptickcl  (displayed  in  the  lower  half  of  the  figure)  but 
to  a  greater  extent  by  tlie  time-averaging  (in  this  case  5  seconds).  Clearly,  an  inverse  relationship  exists 
between  the  degree  of  time-averaging  that  is  employed  and  the  radial  resolution.  The  azimuthal  resolution 
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Figure  2.2  Image  of  a  synthedc  point  source  loeatcd  .35  km  fmm  the  array  at  a  hack 
a/imuth  of  225°.  Wavenumber  Integration  was  used  to  generate  synthetic  seismograms 
individually  for  each  of  the  25  vertical  component  sensors  in  the  NORESS  array  (il¬ 
lustrated  is  the  synthetic  computed  for  the  center  station).  In  this,  and  all  subsequent 
images,  contour  values  indicate  amplitudes  in  dB  relative  to  the  largest  value  in  the 
image  aixl  cylindrical  prop;igalion  of  scattered  wavefronts  was  assumed. 


of  the  Synthetic  source  is  very  pexx  primanly  because  the  coherence  of  the  surface  waves,  computed  with 
only  a  small  degree  of  numerical  noise,  is  very  high  and  becaase  we  are  employing  envelopes  of  the 
bc'arns 

Energy  is  aliased  away  from  tlx’  .tciiial  location  of  tlie  scatterer  to  Uxrations  which  share  the 
same  delav  time.  By  manipulation  of  equation  2.1.  it  can  be  shown  that  scatterers  which  share  a  common 
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delay  time  r  lie  on  a  curve  described  by: 


T 

Ip,  I  - 


(2.4) 


When  p,  is  greater  than  p,  (eg.  P  to  scattering)  this  describes  an  ellipse  with  one  focus  at  the 
center  of  the  array,  major  axis  pointing  to  the  primary  source  and  eccentricity  proportional  to  p,.  In  this 
synthetic  calculation  the  primary  source  was  located  directly  beneath  the  array  (Pi  =  0)  and  thus  the 
curves  of  constant  r  have  degenerated  into  circles  centered  on  the  array. 

To  test  the  dependence  of  the  aximuthal  and  radial  resolution  on  the  range  of  the  source  from  the 
array  an  additional  suite  of  synthetics  has  been  computed  in  which  three  point  sources  were  placed  12.5. 
25.0  and  37.5  km  from  the  array  at  a  back-azimuth  of  225°.  The  image  obtained  using  these  synthetics  is 
di.splayed  in  Figure  2  i.  Not  surprisingly,  the  azimuthal  resolution  is  reduced  as  distance  from  the  array 
to  the  source  increases.  There  is  no  discerrtable  reduction  of  radial  resolution  with  increasing  range  due 
to  dispersion  of  the  surface  wave  packet  since  the  radial  resolution  is  degraded  to  a  greater  degree  by  the 
time-averaging  (in  this  case  5  seconds).  The  azimuthal  resolution  depends  strongly  on  the  aperture  of  the 
seismic  array.  In  Figure  2.4  is  displayed  the  image  of  the  point  source  used  to  create  Figure  2.2.  however 
in  this  case  (he  array  used,  refene '  to  as  MORESS,  is  simply  the  NORESS  array  with  all  dimensions 
doubled  by  affine  transformation.  By  comparison  of  Figures  2.2  and  2.4  it  is  clear  that  the  azimuthal 
rc.solution  has  improved  roughly  by  a  factor  of  2.  The  greater  spacing  between  the  sensors  in  this  new 
array,  however,  means  that  the  spatial  aliasing  is  more  severe  (note  the  increased  size  of  the  side  lobes 
relative  to  those  in  Figure  2.2). 

The  preceding  images  were  computed  using  synthetics  corrupted  only  by  a  small  degree  of 
numerical  noi.se.  In  Figure  2.5  we  test  the  ability  of  the  imaging  algorithm  to  locate  successfully  a 
secondary  .source  when  the  records  contain  a  signiGcant  amount  of  noi.se.  that,  for  some  reason,  cannot 
be  suppressed  by  filtering.  To  generate  this  image  we  have  added  Gaussian  noise,  uncorrelated  in  space 
and  time,  to  the  synthetics  u.scd  to  generate  Figure  2.2.  The  image  has  suffered  a  significant  amount  of 
corruption  but,  despite  high  noise  levels,  evident  in  the  time  series  displayed  in  the  lower  half  of  this 
figure,  the  second;iry  source  is  still  inferred  to  exist  at  the  correct  location. 

The  imaging  algorithm  assumes  that  all  excitation  of  a  .scattcrer  located  at  Rq,  occurs  at 
r  =  Rn,«p,  (the  scattcrer  is  illuminated  by  a  single  arrival).  This  condition  is  met  by  the  synthetics, 
and  as  a  result  it  is  possible  to  inter  the  geographic  location  of  the  synthetic  soua'c,  not  just  the  time 
of  excitation  and  angle  of  approach.  If  this  Is  not  the  case,  and  the  incoming  wavefi-ain  consists  of  a 
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Figure  2.3  Image  of  three  synthetic  point  sources  located  13,  35  and  53  km  from  the 
array  at  a  back  arimuth  of  225°.  Illustrated  is  the  wavenumber  integration  synthetic 
computed  for  tlic  center  station  in  the  array. 


protracted  sequence  of  arrivals,  then  a  single  scatlerer  will  be  illuminated,  and  likely  emit  energy,  for  a 
longer  time. 

In  theory,  this  should  result  in  a  "tail"  of  energy  extending  outward  from  the  actual  location  of 
the  scaiiercr  along  a  radial  line  tscc  Chapter  3).  In  such  cases  it  is  more  appropriate  to  plot  the  image  as 
a  fuiKiKin  of  time,  not  space. 
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Figure  2.4  Image  of  a  synthetic  point  source  located  35  km  from  the  array  at  a  back 
azimuth  of  225°.  Illustrated  is  the  wavenumber  integration  synthetic  computed  for  the 
center  station  in  the  MORESS  array  (an  array  created  by  dilating  the  NORESS  array 
by  a  factor  of  2). 


2.4  Preliminary  analysis  of  recorded  data 


2  4  a  Analysis  of  single  events 


In  this  section  we  apply  the  technique  to  recordings  of  5  nuclear  explosions  and  1  earthquake 
made  by  the  NORESS  arraj  ...  southern  Norway  (.see  Table  2.1).  To  test  the  resolving  power  of  events 
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Figure  2.5  Image  of  a  synthclic  point  source  located  .15  km  from  the  array  at  a  back 
a/imuth  of  225°  computed  after  adding  a  significant  amount  of  Gaussian  white  noise 
to  each  channel.  Illustrated  is  the  wavenumber  integrauon  synthetic  computed  for  the 
center  station  in  the  array. 


coming  fiom  a  single  azimuth  (and  thus  to  gauge  the  improvement  we  can  expect  by  using  a  widely 
disinbuied  suite  of  events)  an  image  has  been  computed  using  a  single  Semipalatinsk  nuclear  explosion 

(S4^51). 

Although  all  Semipalatinsk  explosions  produce  images  that  arc  remarkably  similar,  we  chose 
this  event  because  it  has  been  previously  considered  by  Gupta  el  al ,  (1990a,b)  who  also  applied  beam- 
correcUon.  and  f-k  analysis  in  their  search  for  large  scaltcrers.  This  gives  us  an  opportunity  to  compare 
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Figure  2.6  This  figure  is  included  to  illustrate  the  effect  of  beam  correction.  Pictured 
in  the  upper  box  is  the  NORESS  center  element  (aO)  recording  of  a  Semipalatinsk 
nuclear  explosion  (84351;  see  Thble  2.1).  The  middle  box  displays  the  residual  aO 
channel  after  the  beam  to  the  primary  source  has  been  coherently  subtracted.  The 
residual  channels  have  been  bcamform^  for  a  back  azimuth  of  80°  and  a  phase  velocity 
of  3.29  kms  The  resultant  beam  has  been  converted  to  an  envelope.  Energy  is 
concentrated  at  a  delay  time  r  of  roughly  3  s. 


both  techniques.  Because  of  the  energetic  nature  of  the  Semipalatinsk  tcicseism.  it  has  been  necessary  to 
apply  beam  correction  to  the  recorded  data.  To  illustrate  the  need  for.  and  the  effect  of,  beam  correction, 
we  have  iiKludcd  Figure  2.6.  Illustrated  in  the  upper  part  of  this  figure  is  a  single  channel  recording 
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Figure  2.7  Image  of  the  scaticnng  field  in  the  vicinity  of  the  NORESS  array  obtained 
using  a  single  Semipalatinsk  nuclear  explosion  (843.S1),  The  image  was  generated  by 
employing  1  second  of  time  averaging.  The  test  site  is  IcKaled  at  a  back  azimuth  of 
75°.  Due  to  the  protracted  form  of  tlic  incoming  wavclicld  (estimated  by  the  beam 
pictured  in  the  lower  half  of  the  figure),  the  image  has  been  plotted  as  a  function  time. 
In  practice  we  found  that  applying  a  gain  to  correct  for  geometrical  spreading  of  the 
recorded  scattered  waves  resulted  in  an  unacceptable  amount  of  noise  amplification  at 
extreme  ranges  so  we  did  not  do  it. 


(aO)  of  the  Scmipalatin.sk  nuctciir  explosion  (84.151).  In  the  middle  Ifiird  of  the  figure  is  the  residual 
rcct'rd  that  remains  after  correcting  (coherently  subtracting)  primary  .source  beam  source  from  channel 
;i0.  To  obtain  the  lower  third  of  the  figure  all  array  records  have  been  beam -corrected,  beamformed 
for  a  back-azimuth  of  80"  (an  azimuth  known  to  contain  an  apparent  secondary  source  of  energy).  The 


resultant  coherent  beam  has  been  converted  to  an  envelope.  There  Ls  clearly  a  concentration  of  energy 
centered  at  r  s  behind  the  compressional  onset. 


Thblc  2.1.  Event  Locations 

Event 

Latitude 

Longitude 

Origin  Time 

(°N) 

CE) 

year/day/h;m:s 

84351 

49.88 

78.82 

84/351/03:55:02.8 

85134 

-10.72 

41.26 

85/134/13:25:01.2 

86318 

37.10 

-116.05 

86/318/16:00:00.0 

87225 

37.06 

-116.05 

87/225/14:00:00.1 

87267 

37.23 

-116.38 

87/267/15:00:00.0 

87287 

37.09 

-116.05 

87/287/14:00:00.1 

The  steps  illustrated  in  Figure  2.6  are  systematically  followed  for  all  1681  sub-areas  to  construct 
an  image  using  recorded  data.  The  image  in  Figure  2.6,  and  the  following  ones  computed  using  recorded 
data,  have  been  plotted  as  functions  of  time  because  of  the  protracted  nature  of  the  incoming  primary 
wavetrains.  The  first  image,  obtained  from  the  Semipalatinsk  event  (84351).  is  displayed  in  Figure  2.7  and 
has  been  computed  using  1  second  of  time  averaging.  It  contains  significant  amounts  of  energy  smeared 
along  ellipses  of  constant  r.  Not  surprisingly,  the  ability  of  a  single  primary  event  to  resolve  structure  in 
the  azimuthal  direction  is  poor.  As  in  synthetic  tests,  the  smearing  is  due  to  the  high  coherence  of  the  low 
frequency  scattered  phases  across  the  array.  We  can  decrease  the  smearing  by  increasing  time-averaging 
so  that  several  cycles  of  less  coherent  scattered  energy  are  included  in  the  window.  By  increasing  the 
time  span  to  5  seconds  of  coda  centered  on  the  predicted  onset  we  obtain  the  result  displayed  in  Figure 
2.8.  The  smearing  Ls  still  present,  but  has  been  significantly  reduced,  and  two  apparent  energy  sources 
are  clearly  resolved  although  with  somewhat  lesser  radial  resolution.  The  feature  in  Figure  2.8.  close 
to  the  array  in  the  first  quadrant  is  due  to  the  energy  displayed  in  the  lower  third  of  Figure  2.6  at  r 
s  behind  the  primary  onset.  Previous  studies  (Gupta  e(  ai.  1990a;  Bannister  et  al .  1990)  which  have 
used  the  NORESS  array  to  locate  local  sources  of  scattered  energy  have  also  identified  an  apparent 
local  secondary  source  to  the  south-west  of  the  array.  These  studies  have  concluded  that  this  energy  is 
most  likely  due  to  scattering  interactions  occurring  at  the  Skreikampen/Lake  Mjosa  topographic  relief 
(a  north-west  to  south-east  trending  lake  with  an  adjacent  mountain  that  has  roughly  KXX)  m  of  total 
relicO  located  roughly  35  km  from  the  array.  In  addition,  the  latter  study  noticed  energy  coming  from 
the  north-east  of  the  array  and  attributed  this  to  the  Bronkeberget  topographic  relief  at  a  range  of  roughly 
10  km. 
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To  expedite  the  present  analysis,  we  assume  that  two  point  sources  are  present  and.  as  is 
suggested  by  Figure  2.8  and  the  two  previous  studies  (Gupta  et  at.,  1990a;  Bannister  et  at.,  1990),  are 
located  roughly  3 S  km  to  the  southwest  of  the  array,  and  10  km  distant  at  a  back  azimuth  of  80°.  Since 
we  have  fixed  the  secondary  source  locations,  it  is  possible  to  use  equation  2.4  to  predict  where  energy 
originating  at  these  locations  should  be  aliased  to  in  this  image.  We  have  done  so  and  have  overlain  the 
ellipses  associated  with  these  two  apparent  sources  on  the  image  in  Figure  2.8.  The  strong  coincidence 
between  these  ellipses  and  locally  high  power  levels  suggests  that  much  of  the  energy  in  this  figure  likely 
originates  from  the  two  sources. 

In  Figure  2.9  we  take  this  aitalysis  one  step  further  by  placing  two  point  sources  at  these 
locations,  illuminating  them  using  arrival  times  consistent  with  a  telcseism  arriving  from  Semipalatinsk, 
and  computing  synthetics  for  each  element  in  the  array. 

Although  this  has  been  a  crude  forward  modelling  exercise  (the  synlhctic  the  sources  are 
infinitely  small,  perfectly  omnidirectional  and  impulsively  excited)  this  image  bears  a  striking  resemblance 
to  the  real  image  in  Figure  2.8  despite  the  use  of  a  gain  to  compensate  for  geometrical  spreading  only  in 
the  synthetic  image.  In  the  synthetic  image,  since  only  two  point  sources  exist,  and  their  locations  are 
known,  wc  can  disenminate  between  the  true  sources,  and  false  sources  that  arc  due  to  spatial  aliasing. 
For  example,  considering  the  inner  ellipse,  wc  know  that  the  true  source  exists  to  the  east  or  the  array 
and  that  all  souaes  inferred  to  exist  close  to  the  array  to  the  south,  west  and  north  are  false,  and  are 
onl>  due  to  energy  aliased  away  faim  the  true  location.  The  two  small  sources  inferred  to  exist  in  the 
Norwegian  crust  to  the  south  and  west  of  the  array  (at  a  time  range  of  roughly  3  s)  by  the  image  in 
Figure  2.8  arc  almo- 1  certainly  only  due  to  aliasing  from  a  single  source  to  the  east  of  the  array. 


2  4  h  Image  instability  -  should  the  fold  of  coverage  be  increased’’ 

Clearly  it  would  he  desirable  to  produce  more  robust  images  of  the  l(x;al  scaltcrcrs  by  simulta¬ 
neously  prtKCssing  a  broadly  distributed  suite  of  primary  events.  In  doing  so.  wc  cinild,  in  effect,  image 
the  crust  through  a  "synthetic  aperture  array".  Such  an  event  suite  would  provide  a  more  “balanced" 
illumination  of  the  population  of  scaltcrcrs  that  could  be  provided  by  a  single  event  since  wc  could  scan 
the  crustal  volume  for  stable  apparent  sources.  As  dtscussed  earlier,  when  considering  event  84351, 
it  was  noticed  that  the  azimuthal  resolution  was  improved  somewhat  at  the  expense  of  reducing  radial 
resolution  -  or  in  other  words  -  that  a  trade -off  exists  between  azimuthal  and  radial  resolution  that  cannot 
be  avoided  when  using  single  primary  events.  It  is  possible  to  increase  azimuthal  re.solution  using  a 


40 


-15  -10  -5  0  5  10  15 

time  east  from  NORESS  (s) 


Figure  2.8  Image  of  the  scattering  field  in  the  vicinity  of  the  NORESS  array  obtained 
using  a  single  Semipalatinsk  nuclear  explosion  (84351).  The  image  was  generated  by 
employing  5  seconds  of  time  averaging.  The  beam  computed  for  the  test  site  is  pictured 
in  the  lower  half  of  the  figure.  Overlain  on  the  image  arc  two  ellipses  representing 
the  curves  of  constant  delay  time  r  appropriate  for  scatterers  illuminated  by  this  event. 
The  secondary  sources  are  located  38  km  from  NORESS  at  a  back  azimuth  of  225° 
and  13  km  at  a  back  azimuth  of  80°  . 


broadly  distributed  suite  of  primary  events  without  degrading  radial  resolution.  It  is  certain  that  stacking 
numerous  events  will  help  suppress  incoherent  noise  in  the  event  records,  and  thus,  for  this  reason,  image 
stability  should  increase.  The  question  we  want  to  address  next  is  whether  it  is  feasible  to  attempt  the 
stacking  of  dissimilar  primary  events  to  find  a  more  robust  image  of  the  local  scatterers.  In  an  attempt 
to  answer  this  question,  we  have  computed  an  image  of  the  scatterers  by  using  telcseisms  horn  3  nuclear 
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Figxire  2.9  Image  of  two  synthetic  point  sources  (located  35  km  from  the  array  at  a 
back  azimuth  of  225°  and  13  km  from  the  array  at  a  back  azimuth  of  80°).  The  sources 
were  illuminated  by  an  impulsive  synthetic  primary  event  located  near  Semipalatinsk. 
Illustrated  is  the  wavenumber  integration  synthetic  computed  for  the  center  station  in 
the  array. 


explosions  detonated  at  the  Nevada  Test  Site  (NTS).  The  result  is  displayed  in  Figure  2.10  with  the 
smearing  ellipses  computed  for  the  two  local  sources  considered  in  Figure  2.9  overlain.  The  ellipses 
have  changed  since  p,  has  changed  (the  back  azimuth  to  the  NTS  is  roughly  318°).  A  secondary  source 
is  still  inferred  to  exist  to  the  east  of  the  array  (at  a  time  offset  of  roughly  4  s),  but  beyond  that  coin¬ 
cidence.  this  image  is  profoundly  different  from  that  in  Figure  2.8.  A  strong  source  Ls  inferred  to  exist 
to  the  northwest  of  the  array.  Several  coherent  sources  are  suggested  to  exist  to  the  south  of  the  array. 
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Figure  2.10  Image  of  the  scattering  field  in  the  vicinity  of  the  NORESS  array  obtained 
using  four  nuclear  explosions  detonated  at  the  Nevada  Ibst  Site  (NTS)  (86318.  87225, 
87267  and  87287).  The  image  was  generated  by  employing  5  seconds  of  time  averaging. 
The  beam  computed  for  the  test  site,  using  recordings  of  the  second  event,  is  pictured 
in  the  lower  half  of  the  figure.  Overlain  on  the  image  are  ellipses  representing  the 
curves  of  constant  delay  time  r  appropriate  for  scatterers  illuminated  by  these  events. 
The  secondary  sources  are  located  38  km  from  NORESS  at  a  back  azimuth  of  225° 
and  13  km  at  a  back  azimuth  of  80°  . 


These  two  images  suggest  that  there  is  the  possibility  that  the  intrinsic  appearance  of  the  crust  depends 
on  the  direction  from  which  it  is  illuminated.  If  we  accept  the  hypothesis  that  topographic  features  are 
responsible  for  the  bulk  of  the  scattered  energy  then  an  examination  of  this  surface  quickly  reveals  that 
common  topographic  structure  has  geographic  extent  and  is  clearly  not  cylindrically  symmetric  -  there  is 
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Figure  2.11  Image  of  the  scattering  field  in  the  vicinity  of  the  NORESS  array  obtained 
using  recordings  of  a  single  earthquake  (85134;  10.7°  S,  41.3°  E).  The  image  was 
generated  by  employing  5  seconds  of  time  averaging.  The  beam  computed  for  the 
earthquake  is  pictured  in  the  lower  half  of  the  figure.  Overlain  on  the  image  are  ellipses 
representing  the  curves  of  constant  delay  time  r  appropriate  for  scatterers  illuminated 
by  this  event.  The  secondary  sources  arc  kxated  38  km  from  NORESS  at  a  back 
azimuth  of  225°  and  13  km  at  a  back  azimuth  of  80°  . 
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a  clear  northwest-southeast  fabric  in  which  typical  features  ae  lO’s  of  km  long.  In  view  of  the  geographic 
extent,  and  shape,  of  typical  topographic  features  it  is  not  be  surprising  that  there  appears  to  be  some 
dependence  of  the  appearance  of  the  crust  on  the  direction  from  which  it  is  illuminated.  If  this  was  the 
only  possible  explanation  of  the  dissimilar  nature  of  the  images  then  perhaps  a  strong  argument  could  be 
made  against  stacking. 
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Fortunately,  other  explanations  exist  which  can  explain  the  disparate  nature  of  these  images 
and  strongly  support  the  simultaneous  use  of  many  primary  events  to  attempt  a  more  robust  imaging 
analysis.  The  image  quality  of  Figure  2.10  is  clearly  de^.^ded  since  the  signal  to  noise  ratio  (SKR)  is 
much  lower  in  the  NTS  seismograms  than  it  is  in  the  Semipalatinsk  events.  We  expect  that  noise  in  the 
data  will  map  directly  into  image  noise,  or  instability.  Some,  or  all.  of  the  apparent  sources  in  this  image 
(or  any  image  produced  using  a  single  event)  could  be  due  to  coherent  bundles  of  noise.  Furthermore, 
a  distinct  possibility  exists  that  the  energy  kKatcd  in  the  fourth  quadrant  of  Figure  2.10  (at  roughly  the 
back  azimuth  of  the  primary  event)  docs  not  have  a  local  source  but.  in  fact,  is  due  to  incomplete  beam- 
correction.  This  argument  could  be  used  to  shed  doubt  on  the  reality  of  the  secondary  source  inferred  to 
exist  to  the  cast  of  the  array  by  the  Semipalatin.sk  event  (Figure  2.8)  since  Semipalatinsk  is  located  at  a 
back -azimuth  of  75°  from  the  array. 

Obviously,  neither  of  the  two  images  considered  thus  far  can  confidently  be  touted  as  repmsen- 
lativc  of  the  true  distribution  of  scaltcrers.  Unle.ss  the  intrinsic  appearance  of  the  crust  is  wildly  dependent 
on  the  direction  of  illumination  it  seems  clear  that  image  quality  can  be  improved  if  a  large  number  of 
events  are  allowed  to  “vote"  on  what  a  representative  image  of  local  scaltcrers  should  be.  It  also  seems 
logical  that  greater  attention  should  be  paid  to  the  events  recorded  with  higher  SNR.  If  the  feature  in 
the  fourth  quadrant  in  Figure  2.10  is  due  to  incomplete  beam  correction,  and  even  if  this  problem  is 
shared  by  many  events,  then  stacking  of  a  bro;idly  distributed  suite  should  suppress  this  effect,  since  the 
primary  source  energy  should  lie  in  a  pic  slice  originating  at  the  array  and  pointing  toward  each  primary 
source.  Furthermore,  such  slacking  should  serve  to  suppress  the  influence  of  near  source  scatterers  -  and 
incoherent  noise. 

A  very  similar  problem  has  plagued  the  reflection  seismic  community.  As  in  our  work,  they 
image  subsurface  structure  using  recordings  of  phases  produced  by  the  interaction  of  man-made  seismic 
energy  and  sub-surface  reflectors.  It  is  well  known  (Yilmaz.  1987:  Claerbout,  1985)  that  attempting 
such  imaging  using  data  that  has  .sampled  each  point  on  a  reflector  of  interest  only  once,  “single  fold 
coverage”,  generally  yields  poor  results.  Given  the  .size  of  typical  primary  sources,  the  SNR  isn’t  high 
enough  to  allow  good  resolution  of  the  subsurface  “secondary”  sources.  To  improve  the  results,  seismic 
data  is  generally  collected  such  that  redundant  excitation  of  each  subsurface  point  of  interest  is  achieved 
-  or  to  use  the  parlance  of  the  reflection  seismology  community  -  the  fold  of  coverage  is  increased.  The 
redundant  recordings  are  stacked  to  interfere  coherently  such  energy.  In  the  current  problem,  we  can 
increase  the  fold  of  coverage,  or  we  can  in  effect  cause  each  active  scaticrer  to  be  redundantly  excited, 
by  using  more  than  one  primary  event. 
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2  A  c  An  illustration  of  why  temporal  deconvolution  is  necessary 

As  mentioned  in  section  2.2,  further  complication  is  provided  by  fact  that  the  excitation 
futKUon  is  clearly  not  impulsive  but  consists  of  a  protracted  sequence  of  arrivals  (a  wavetrain).  In  Figure 
2.12  is  illustrated  an  image  computed  using  a  single  telcseismic  earthquake  as  a  primary  source  (event 
851.14  located  at  a  back  azimuth  of  133'')  and  a.ssuming  a  pj  of  ..304  skm  '.  An  energy  low  in  the 
second  quadrant  (at  the  back  azimuth  to  the  primary  event)  suggests  that  the  beam  correction  has  been 
effective.  The  two  sources  discussed  earlier  are  still  prominent  but  an  additional  apparent  source  is 
inferred  to  exist  in  the  third  quadrant.  When  generating  images  using  temporally  protracted  teleseisms,  a 
fundamental  ambiguity  exists  when  trying  to  locate  sources  of  scattered  energy.  Is  one  apparent  source 
inferred  to  exist  at  a  greater  distarKC  from  the  array  than  another  because  it  is  located  farther  from  the 
array  or  because  it  is  due  to  excitation  by  a  tardy  primary  arrival?  It  .seems  likely  that  this  distant  source 
IS  simply  due  to  a  late  arrival  from  the  earthquake  exciting  a  scaiterer  located  relatively  close  to  the 
array  since  other  events  sec  no  source  at  this  location.  Prior  to  generating  spatial  images  of  the  local 
scattering  field  using  recorded  telcseismic  events  as  primary  sources  it  is  clear  that  an  attempt  should  be 
made  to  suppress  the  intiuence  of  late  primary  arrivals.  In  Chapter  3  a  simple  convolutional  model  is 
dcvelo(x.'d  to  describe  the  interaction  between  protracted  tele.seismic  wavetrains  and  the  impulse  response 
of  the  liK'al  scatterers.  A  remedy,  in  further  prc-proce.ssing  involving  deconvolution  of  an  estimate  of  the 
inciKTiing  wavetrain  from  the  individual  records,  is  developed. 


2.5  Conclusion 

In  this  chapter,  a  technique  to  image  for  local  sources  of  scattered  energy  using  small  aperture 
iirray  recordings  has  been  developed.  In  essence,  a  hyperboloid  summation  migration  is  used  to  resolve 
iho  faint  secondary  sources  of  R,  energy.  This  technique  tws  been  applied  to  array  recordings  of  individual 
real  and  synthetic  events.  The  synthetics  have  shown  that  the  azimuthal  spatial  resolution  provided  by  the 
technique,  when  applied  to  a  small-aperture  array,  is  inversely  proportional  to  the  distance  to  the  secondary 
source  and  directly  proportional  to  the  aperture  of  the  array.  In  theory,  the  radial  resolution  should  be 
inversely  proportional  to  the  range  to  the  .second;uy  .source  (due  to  dispersion  of  the  surface  wave  p  icket). 
Hossever  this  effect  is,  in  practice,  lessened  due  to  the  effect  of  time-averaging.  Preliminary  analysis  of 
several  events  recorded  by  the  NORESS  array  in  southern  Norway  indicates  that  stable  secondary  sources 
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likely  exist  in  the  vicinity  of  the  array  (to  the  southwest  and  northeast).  These  results  are  consistent  with 
those  of  Gupta  et  a/  (i990a,b)  and  Bannister  et  <j/.(l990).  It  appears,  however,  that  the  production  of 
robust  images  cannot  likely  be  achieved  using  single  primary  events.  The  feeble  secondary  sources 
cannot  be  well  imaged  without  redundant  excitation.  It  is  necessary  to  “increase  the  fold"  by  stacking 
the  results  of  a  large  suite  of  broadly  distributed  primary  events.  Because  of  the  protracted  nature  of 
most  teleseisms  it  is  not  clear  what  primary  phases  are  responsible  for  exciting  locally  scattered  energy. 
As  a  result,  only  temporal  images  have  been  produced.  To  refine  the  Jgorithm  further,  to  facilitate  the 
production  of  spatial  images,  will  require  deconvolution  of  the  incident  wavetrain  (Chapter  3).  This  will 
allow  us  to  test  more  stringently  the  stability  of  our  images. 
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Chapter  3 


Regularized  deconvolution  of  noisy  seismic  array  recordings 

Abstract 

In  this  chapter,  an  algorithm  is  develt^d  for  the  deconvolution  of  a  presumed  known  wavelet 
from  uncertain  data.  This  approach  seeks  a  deconvolved  trace  (model)  that  possesses  a  physically 
reasonable  amount  of  covariance  and  produces  an  acceptable  data  misfit.  Regularization  is  achieved 
using  the  L2  norm  to  assess  data  misfit  and  model  roughness  (or  complexity).  The  method  of  Lagrange 
multipliers  is  used  to  control  the  tradeoff  between  these  2  characteristics.  The  noise  is  assumed  to  result 
from  a  stationary,  Gaussian,  process.  The  noise  covariance  is  estimated  assuming  the  process  is  ergodic 
and  applying  an  autocorrelation  operator  to  pre-event  noise. 

Two  approaches  are  taken  to  introduce  the  proper  level  of  covariaiKe  into  the  model.  In  one, 
akin  to  the  Occam’s  Razor  method  of  Constable  et  a/.(l987).  a  straightforward  differencing  operator  is 
applied  to  the  model  vector  to  yield  a  smoothed  result.  In  the  second  approach,  a  priori  knowledge, 
obtained  from  physical  arguments,  is  used  to  provide  an  estimate  of  the  statistical  nature  of  the  model. 
This  information  is  used  to  define  a  correlation  operator  that  is  applied  to  the  model  vector  as  a  way  of 
seeking  a  model  with  an  appropriate  degree  of  covariance.  In  both  approaches,  it  is  assumed  that  the 
mailing  kernel  which  relates  the  model  to  the  observed  data  is  known  exactly.  A  series  of  synthetic 
models  is  explored  in  an  attempt  to  define  weaknesses  of  the  approaches.  Comparison  is  made  with  the 
spectral  division  technique. 

Full  array  synthetics  are  convolved  with  recorded  time  scries  and  then  deconvolved,  after  adding 
nui.se,  to  determine  the  extent  to  which  the  resolution  of  cru.stal  .scaltercrs  neat  seismic  arrays  c;m  be 
cnhanecd  by  using  this  processing. 
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C.l  Introduction 


The  imaging  analysis  introduced  in  Chapter  2  is  based  on  equation  2.3 


r  =  Ra. -(p, -pj  (3.1) 

that  assumes  impulsive  excitation  of  the  local  population  of  scatterers  to  predict  arrival  times  of  scattered 
phases  at  seismic  receivers.  This  formula  was  applied  to  recordings  of  teleseismic  events  to  produce 
images  of  the  scattering  field  in  the  vicinity  of  the  NORESS  array.  If  one  adopts  the  simple  assumption 
that  the  level  of  excitation  of  an  active  scattercr  is  proportional  to  the  amplitude  of  the  incoming  wavefield 
then  it  is  clear  that  excitation  is  not  impulsive.  Instead,  due  to  the  protracted  nature  of  the  incoming 
wavefields,  estimated  by  beamforming,  scatterers  are  repeatedly  illuminated  by  successive  arrivals  in  a 
typical  teleseismic  wavetrain.  Given  this  hypothesis,  it  was  concluded  in  Chapter  2  that  formula  2.3 
should  be  used  to  produce  temporal,  not  spatial,  images  of  the  local  scatterers.  It  is  clear  that  further 
pre-processing  of  the  array  records,  with  the  goa.  of  reducing  the  effective  time  span  of  teleseismic 
arrivals,  is  desirable. 


C.2  Analysis  of  the  forward  problem 


C.2.a  The  convolutional  model 


In  an  attempt  to  remedy  this  situation  we  assume  a  simple  model  for  the  array  records.  The 
motions  occurring  at  the  locatirxi  of  an  individual  sensor,  d(l),  are  real  values  that  result  from  the 
convolution  of  the  sensor’s  impulse  response,  m(t),  hereafter  known  as  the  model,  and  the  sequence 
of  arrivals  from  the  primary  source,  b(t).  which  we  assume  is  invariant  with  respect  to  time  shift.  In 
more  precise  terms,  we  define  the  impulse  response,  m(t),  to  be  the  wavefield  that  would  be  recorded 
by  that  sensor  in  the  event  of  impulsive  excitation  by  a  plane  wave  with  the  wavevector  appropriate  for 
the  teleseismic  source.  The  wavetrain,  b(t),  goes  under  many  aliases  -  the  transfer  function,  the  basis 
wavelet,  the  mapping  kernel,  the  representer,  the  data  kernel,  etc.  Hereafter,  in  this  chapter,  b(t)  will 
be  referred  to  simply  as  the  wavelet.  Considering,  for  the  moment,  recordings  uncomipted  by  additive 
noLse: 


d{t)  =  b(l) 


h(t  —  T)m{T)dT  = 


•  i 


*2 

j  g{t,T)m(T)dT 


(3.2) 
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where  *  denotes  convolution.  The  model  is  assumed  to  have  finite  support  between  ti  and  t2.  g(<,  T)=b(t- 
r).  It  is  assumed  that  these  motions  superpose  linearly  (and  that  we  have  a  linear  forward,  and  thus 
linear  inverse,  problem).  Given  that  we  consider  discrete  representations  of  the  data,  model  and  wavelet 
these  relations  can  be  rewritten  as: 

dk  =  =  ^9k}mj  (3.3) 

j-i  j-i 

where  gkj  =  bk-j  +  \  is  defined  from  j=l  to  i=Ni,,  the  number  of  points  in  the  wavelet.  The  number  of 
model  components  is  given  by  the  data  is  recorded  at  points.  Given  this  discrete  representation, 
the  solution  of  the  error-free  linear  forward  problem  can  be  rewritten  as  a  matrix  equation 

d  =  Gm  (3.4) 


where  G  is  a  matrix  that  contains  the  wavelet  b(t)  and  can  be  defined  as  follows: 


G  = 


/  6,  0 

6.  6,  0 


bn  - 1  ...  6)  0 

0  b„  6„_i  ...  b\  0 


0  \ 


0  bn 


be  bi 


(3.5) 


I .  0  bn  bn  I  I 

Vo  .  0  b„  / 

in  which  ca.se  the  mauix  is  not  square  (N^  =  Nm  +  Nb  -  !)•  Given  that  N6  is  greater  than  one.  there 
arc  more  data  points  than  model  parameters  to  be  determined.  The  problem,  as  defined  above,  is  likely 
over  determined.  All  the  direct  constraint  information  we  have  to  form  an  estimate  of  the  model  is 
contained  within  Equation  C.4.  All  wavelets,  g^.  J=l,  Nd  are  linearly  independent  and  thus  form  the 
basis  for  an  Nd  dimensional  linear  vector  space,  known  hereafter  as  the  data  space.  A  subspace  of  the 
data  space  encompassing  the  first  N„,  points  is  hereafter  referred  to  as  the  model  space. 


C  2  b  The  effect  of  protracted  teleseismic  wavelets  on  imaging  for  local  scatterers 

Prior  to  considering  the  problem  of  inverting  equation  C.4.  it  seems  reasonable  to  apply  the 
convolutional  model  to  synthetic  and  recorded  data  and  to  assess  the  effect  of  a  prolonged  excitation 
funcuon  on  the  imaging  process.  In  Figure  C.l  is  displayed  the  image  computed  when  a  single  scattercr 
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Figure  3.1  Image  of  a  synlhetic  point  source  located  35  km  from  the  array  at  a  back 
azimuth  of  225°.  Wavenumber  Integration  was  used  to  generate  synthetic  seismograms 
individually  for  each  of  the  25  vertical  component  sensors  in  the  NORESS  array  (il¬ 
lustrated  is  the  synthetic  computed  for  the  center  station).  In  this,  and  all  subsequent 
images,  contour  values  indicate  amplitudes  in  dB  relative  to  the  largest  value  in  the 
image  and  cylindrical  propagation  of  scattered  wavefronts  was  assumed. 


35  km  to  the  southwest  of  the  NORESS  array  is  excited  by  a  single,  impulsive,  arrival  impinging  on  the 
area  from  directly  beneath  (it  has  infinite  phase  velocity). 

The  radial  resolution  is  degraded  i»nly  by  dispers-jn  and  time-averaging  which,  in  this  case,  is 
5  seconds.  The  azimuthal  resolution  of  the  NORESS  array  at  this  range  is  roughly  20  km,  the  image 
sidelobes  are  not  significant.  In  Figure  C.2.  we  display  a  recorded  lele.seism,  produced  by  a  Semipalatinsk 
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Figure  3.2  A  Scmipalatinsk  tclescismic  event  recorded  by  the  NORESS  array.  Dis¬ 
played  in  (be  upper  figure  is  (he  beam  computed  using  24  channels,  the  lower  figure 
displays  the  amplitude  spcctnim  of  this  beam. 


nuclear  explosion  (1985,  day  351).  The  character  of  the  wavcficld  arriving  from  the  explosion  has  been 
estimated  using  simple  delay  and  sum  beamforming  -  we  denote  this  beam  us  b,,(l).  To  produce  the 
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Figure  3.3  Displayed  is  the  image  computed  using  time  series  that  have  resulted  from 
the  convolution  of  the  synthetics  displayed  in  Figure  C.l  and  the  recorded  teleseism 
displayed  in  Figure  C.2.  Illustrated  in  the  lower  figure  is  the  aO  synthetic/recorded 
seismogram.  This  figure  is  included  to  illustrate  what  degradation  of  radial  resolution 
we  expect  when  locid  point  scatterers  are  excited  by  a  relatively  impulsive  tcleseism. 


image  in  Figure  C.3.  we  have  used  this  teleseism  to  play  the  role  of  the  wavelet  (excitation  function)  - 
the  teleseism.  bj(t),  has  been  convolved  with  the  synthetics  used  to  produce  Figure  C.l. 

In  the  lower  part  of  this  figure,  a  sample  seismogram  (station  aO  -  vertical  component)  is 
displayed  -  the  convolution  has  clearly  increased  the  duration  of  the  arrivals.  As  a  result,  the  image  of 
the  point  secondary  source  (displayed  in  the  upper  part  of  Figure  C.3)  has  acquired  a  “shadow”  extending 
outward  from  the  array  originating  at  the  true  source  of  scattered  energy  -  the  radial  resolution  has  been 
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Figure  3.4  A  tclescismic  earthquake  (1985,  day  134)  recorded  by  the  NORESS  array. 
Displayed  in  the  upper  figure  is  the  beam  computed  using  25  channels,  the  lower  figure 
displays  the  amplitude  spectrum  of  this  beam.  This  event  is  the  same  one  that  appears 
in  Figure  3.11. 
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Figure  3.5  Displayed  is  the  image  computed  using  time  series  that  have  resulted  from 
the  convolution  of  the  synthetics  displayed  in  Figure  C.l  and  the  recorded  teleseism 
displayed  in  Figure  C.4.  Illustrated  in  the  lower  figure  is  the  aO  synthetic/recorded 
seismogram.  This  figure  is  included  to  illustrate  what  degradation  of  radial  resolution 
we  expect  when  local  point  scatterers  are  excited  by  a  relatively  protracted  teleseism. 
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“pulsy"  nature  of  the  earthquake  wavetrain  in  Figure  C.4  has  translated  directly  into  the  image  as  discrete 
concentrations  of  energy  located  mainly  to  the  southwest  of  the  true  scatterer.  Still  plotted  as  a  function 
of  space,  this  image  implies  the  existence  of  a  fai.se  source  roughly  70  km  to  the  southwest  of  the  array. 
We  know  that  this  energy  originated  only  3.^  km  from  the  array  but  is  due  to  a  multiple  that  arrived 
roughly  10  seconds  after  the  primary  onset  The  effect  of  the  telcseismic  nuclear  explosion  on  the  image 
of  the  local  scatterers  seems  relatively  light  and.  given  that  the  convolutional  model  gives  a  reasonable 


description  of  how  the  energy  in  the  wavelet  interacts  with  the  scatterers,  should  be  relatively  easy  to 
remove.  The  earthquake,  however,  is  characterized  by  more  significant  coda  which  will  be  harder  to 
correct.  In  these  synthetic  examples,  the  temporal  resolution  of  the  array  records,  and  thus  the  radial 
resolution  of  the  images,  is  being  degraded  by  the  protracted  nature  of  the  wavelet.  Given  that  the 
convolutional  model  introduced  at  the  beginning  of  this  section  gives  a  reasonable  description  of  the 
process  at  work,  it  seems  likely  that  a  similar  process  degrades  the  quality  of  the  images  presented  in 
Chapter  2.  The  goal  of  this  chapter  is  to  devel<^  a  means  by  which  the  radial  resolution  of  the  images  can 
be  enhanced  without  degrading  the  azimuthal  resolution.  The  quality  of  images  produced  using  recorded 
data  will  be  increased  if  the  basic  seismic  wavelet  can,  in  effect,  be  compressed  into  a  delta  function  so 
we  can  use  the  impulse  response  to  image  the  local  population  of  scatterers.  The  goal  of  the  next  section 
is  to  derive  a  means  whereby  the  inverse  probler  an  be  solved.  Basically,  the  matrix  equation,  C.4, 
has  to  be  inverted,  or  the  wavelet  within  G  is  to  be  deconvolved  from  the  data  vector,  d,  to  yield  an 
estimate  of  the  model,  th. 


C.3  Inversion  methodology 

As  we’ve  seen  above,  the  convolution  of  the  impulse  response  of  a  systent  with  a  non-impulsive 
function  leads  to  a  loss  of  resolution,  or  a  blurring  of  detail.  A  large  number  of  attempts  have  been  made 
to  reverse  this  process,  or  deconvolve  the  data  for  the  underlying  impulse  response.  This  section  describes 
three  approaches. 


C  J.a  Deconvolution  by  spectral  division 

One  technique  that  has  been  used  involves  a  straightforward  use  of  the  convolution  theorem. 
Kcprc,scnling  anguhtf  frequency  by  u.  if  D(u;),  B(u;)  and  M(u;)  arc  Fourier  translorm  pairs  of  d(t),  b(t) 
;ind  m(t)  respectively  then  by  Fourier  transforming  equation  C.2  we  find  that: 

D{u})  =  (3.6) 

An  estimate  of  the  spectrum  of  the  model,  M (w)  can  be  found  using  spectral  division: 

A>(w)  =  /)(u;)/i  '(v)  (3.7) 
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and  the  underlying  model  estimate,  rh(l)  can  be  found  by  inverse  Fourier  transformation.  At  first  sight, 
this  technique  seems  like  it  should  work  well,  however  in  practice  spectral  holes  exist  in  B(u)  making 
its  inversion  unstable  -  straight  division  amplifies  frequencies  absent  from  B(u)  which  produces  ringing 
in  M(u).  In  practice,  the  inversion  is  stabilized  by  filling  in  the  spectral  holes  (or,  so  to  speak,  raising 
the  water  level)  present  in  B(cli)  prior  to  division  in  a  process  known  as  pre-whitening  (YiJmaz,  1987). 
The  degree  to  which  the  holes  arc  filled  up  is  most  often  determined  by  visual  inspection  of  m{t)  and 
thus  the  technique  is,  in  general,  highly  subjective.  This  approach,  however,  reveals  an  important  aspect 
of  the  inversion  that,  given  a  finite  number  of  inexact  observations,  cannot  be  avoided.  It  is  well  known 
(Parker,  1977)  that  the  problem  of  inverting  a  linear  forward  problem  for  a  discrete,  or  continuous,  model 
given  a  finite  set  of  inexact  observations  is  non-unique.  The  challenge  lies  in  finding  an  inversion  scheme 
that  avoids  the  ad  hoc  selection  of  a  single  model. 


CJ.b  Wiener  filtering  ■  seeking  the  matimum  resolution  and  the  best  daia-misfit 

As  discussed  by  Oldenburg  (1981)  a  powerful  alternative  solution  of  this  linear  forward  problem 
lies  in  the  inversion  formalism  of  Backus  and  Gilbert  (1967,1968.1973).  At  the  heart  of  this  approach 
is  the  lesson  that,  given  the  finite,  inexact,  nature  of  the  observations,  a  single,  exact,  model  cannot  be 
produced.  It  is  possible  only  to  estimate  blurred  versions  of  the  true  model  using  linear  combinations  of 
the  observations.  In  matrix  form  the  latter  part  of  this  statement  can  be  rewritten  as: 


th  =  Ad 


(3.8) 


or  explicitly  defining  the  terms  involved: 

combining  equations  C.3  and  C.8  we  find  that; 


.Vj 

rhj  =- 

I-  I 


letting 


(3.9) 


(3.10) 


(3.11) 
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where  R  is  generally  known  as  the  resolution  kernel,  we  find  that: 


thj  —  ^ (3.12) 

The  resolution  kernel  describes  the  degree  to  which  the  true  nuxlel  is  blurred.  It  could  be  argued  that  we 
might  want  to  seek  those  data  weights,  a,  which  yield  the  greatc.sl  resolution,  or  the  most  time  limited 
resolution  kernel.  The  greatest  resolution  that  the  data  allow  can  be  found  by  making  a  as  close  as 
possible,  in  the  least  squares  sense,  to  the  delta  function  This  result  can  be  found  by  minimizing; 

N...  Wd 

t-1  k~\ 

with  respect  to  the  data  weights,  a.  It  can  be  shown  that,  provided  G  is  non-singular,  this  yields  the 
following  solution,  expressed  in  matrix  form; 


A  -  (G'G)  'G'’  (.3.11) 

This  set  of  weights  al.so  yields  the  mtxlel  that  has  the  smallest  Euclidean  norm,  or  is  smallest  in  the  least 
squares  sense  and  yields  the  smallest  data  misfit,  as  measured  by  the  method  of  least-squares  -  it  is  a 
minimum  eiror  v:uuuKe  technique  (Mendel,  198.31.  This  result  is  known  as  the  unconstrained  Wiener 
shtqhng  liher  i  Ireitcl  and  Robinson,  1966;  Oldenburg.  1981  and  Mcnkc.  1984).  In  the  unlikely  event 
that  the  observations  are  uncomipted  by  noise  the  predicted  data,  d^,  arc  given  by 

d,,  -  G(G'G)  ’G'<1  (;J.1,5) 

which,  by  inspection,  shows  that  the  data  arc  predicted  exactly.  The  matrix  which  prcmultiplics  the 
observed  data  vector  is  known  as  the  data  resolution  matrix  (Ri;  Mcnkc.  1984)  and,  in  this  case,  equals 
I,  the  identity  matrix. 


(  '  (  An  approach  to  ila  iinvoluiion  that  allows  for  a  ttadcoff  between  data  misfit  and  model  simplicity 


Them  are  several  reasons  why  the  cunent  imaging  analysis  should  benefit  from  a  less  specific 
approach  In  practice,  it  is  always  true  that  the  observations  lue  corrupted  by  noise.  An  equation  that 
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gives  a  more  rc'ilisUc  description  of  the  forward  problem,  than  is  provided  by  equation  C.4.  is 

d  —  Cm  n  (3  16) 

where  the  vector,  n,  rvprc.sents  the  additive  noise.  We  need  an  approach  that  will  yield  a  model  that 
produces  the  correct  level  of  data-misfit  (given  data  made  inexact  by  noise,  an  exact  fit  is  no  longer 
desirable).  Furthermore,  it  is  appropriate  to  cast  the  solution  in  a  way  that  will  allow  the  introduction  of 
a  priori  information,  such  as  that  which  describes  the  expected,  or  desired,  character  of  the  model. 

As  in  many  other  inverse  problems,  the  model  we  seek  suffers  from  a  number  of  undesirable 
qualities  that  we  want  to  suppress  simultaneously  (in  this  case  model  complexity  and  the  inability  to 
reproduce  the  input  data).  A  constrained  optimization  that  will  allow  us  to  control  model  complexity  and 
vary  the  fidelity  with  which  we  reproduce  the  input  data  can  be  achieved  using  the  method  of  Lagrange 
multipliers.  Following  this  approach,  we  define  the  trade-off  functional,  U,  as  follows: 

U  =  ||Om||2  +  /i  ■ '  {IlGm  -  d||2  -  X^}  (3.17) 

where  the  norm  used  is  the  standard  Euclidian  norm,  O  is  an.  as  yet.  undefined  operator  that  wc  apply 
to  the  model  vector  and  X  ^  is  the  mean  square  misfit.  One  reason  why  this  approach  is  attractive  is  that 
it  is  general,  in  the  senc^  ^ve  arc  free  to  use  a  priori  information  to  define  the  norms  in  the  way  we 
feel  best  suits  the  problem  at  hand. 

Considering  the  problem  which  has  motivated  the  work  in  this  chapter,  there  are  several  sources 
of  noise.  In  addition  to  cultural  sources,  local  wind,  and  spring  runoff,  an  important  source  of  noise 
in  the  vicinity  of  the  NORESS  array  is  microseismic  noise,  generated  by  interactions  between  swell 
propagating  in  opposite  directions,  for  example  in  the  North  Sea  (Fyen,  1986).  The  wind  energy  is 
converted  to  acoustic  energy  within  the  water  column  and  then  to  seismic  energy  at  the  seafloor.  The 
dominant  period  of  microseismic  noise  .shows  some  variation  (Hedlin  and  Orcutt,  1989)  but  tends  to 
lie  in  the  band  from  3  to  10  seconds.  Given  that  the  NORESS  array  rec  Is  are  digitize^,  at  40  Hz, 
it  is  immediately  apparent  that  the  noise  in  this  dataset  must  Iv'  correlated  and  the  inversion  theory 
should  accommodate  this  complexity.  Considering  the  complex  nature  of  the  noise,  that  many  sources 
arc  simultaneously  at  wtxk,  wc  assume  the  added  vector,  n.  is  a  single  realization  of  a  stationary  random 
procc.ss.  N(l),  that  obeys  Gaussian  statistics.  By  dclinition,  assuming  N(t)  has  been  demeaned,  the  noise 
autocovariance  is  defined  as; 

/f(T)  =  t-{,V(0.A(t-r)} 
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(3.18) 


where  E  denotes  the  expectation  operator.  Since  we  have  only  a  single  realization  of  N(t)  in  the  pre-event 
noise  sample  we  assume  ergodicity  and  estimate  R(r)  by  a  simple  delay  and  sum  operation  over  J  points: 


1 

J 


J-j 

^7i(t,)n(t,  -t-  jAt) 
1-1 


(3.19) 


and  use  this  infwmation  to  construct  a  Tocplitz  aulcKovariance  matrix,  C,„.,  where  the  (i,j)‘^  element  is 
given  by  R{{i  -  j)At).  We  assume,  for  the  time  being,  that  no  error  exists  in  G,  and  thus,  give,  the 
assumed  Gaussian  nature  of  N(t),  that  an  appropriate  measure  of  data  misfit  u.ses  the  Euclidean,  or  L2. 
norm; 

llGm  -  dl!^  =  (Gm  -  d)’’  (Gm  -  d)  (3.20) 


This  misfit  measure,  known  as  the  method  of  weighted  least  squares  (Priestley,  1981),  is  essentially  the 
one  used  by  Constable  et  uf  (1987).  However,  in  the  present  work  the  noise  is  allowed  to  possess  an 
arbiUury  degree  of  covariance.  .As  argued  by  Constable  ef  a/.(1987),  in  the  event  the  errors  result  from 
a  zero-mean,  Gaussian,  prcxress  and  is  uncorrelated  the  expected  misfit,  after  normalization,  is  1.  As 
illustrated  by  Joumcl  (1989).  this  is  still  true  when  the  errors  are  correlated.  The  covariance  matrix  Cr.n. 
operating  in  this  manner  within  the  norm,  rotates  the  correlated  vector  Gm  -  d  such  that  the  norm  can 
be  rewritten  as  the  dot  product  of  a  new  uncorrelated  vector  which  has  a  covariance  matrix  that  equals 
I,  the  identity  riiainx.  Thus  the  expected  rms  misfit  is  still  1.  The  statistics  underlying  the  model  are  not 
so  easily  aiciis v.ed  -  when  dealing  with  recorded  data,  no  model  sample  exists  (it  is  what  we  seek)  and 
so  m(x]el  Statistics  must  be  inferred  using  Ic.ss  direct  means.  In  one  approach  that  has  been  taken  to  this 
type  of  problem  (i  e.,  Shaw  and  Orcutt,  iqs.*!;  Constable  et  al .  1987),  a  simple  (smooth)  model  is  sought 
by  .ipplymg  a  “roughening  .•nalrix",  R,  to  the  mixlel  vector.  Tlic  roughening  matrix  is  generally  a  first, 
or  second,  order  discrete  differencing  operator.  Following  these  authors,  we  de  inc  the  Lt  model  norm 
as  tiillows. 

I  iii|  =  (Rm)^  (Rni)  —  m^  Rm  (3.21) 


Expanding  the  right  side  of  equation  C.17.  after  sebslituling  equations  C.20  and  C.21: 


ii/H^Hiii-//  '(ni'G'C„,;Gni  2m '  G  '  C’,„;  .1  .  d  '  C,,,' d  A 


1  J  /  I  ^  f  '  J  I  ^  ^ 


(3.22) 


Follow  inp  the  melhtxl  of  l.agnmgc  multipliers,  to  hnd  the  mixlel  vector  that  yields  the  smallest  functional. 


U,  we  differentiate  the  functional  with  respect  to  m  and  force  the  result  to  vanish. 

am 

GC„^d  =  iiR^Rm  +  G^C„^Gm 
(/.R^R  +  G^C„,;g)  m  =  G^C„;.d 
These  steps  yield  a  model  estimate; 

rh  =.  (r^R  +  //  'G^C„^g)  '  /i  ’G^Cjd  (3.23) 

Hereafter  we  refer  to  this  model  estimation  proc''dure  as  method  B.  TTiis  approach  is  attractive  because 
it  guarantees  that  no  unnecessary  structure  is  given  to  u.  model  -  and  thus  it  has  been  named  the 
“Occam’s  Razor”  approach.  The  degree  to  which  structure  is  suppressed  is  determined  by  the  order 
of  the  derivative.  One  limitation  of  this  approach  exists  because  it  does  not  allow  a  rigorous  selection 
of  a  model  roughening  operator.  It  simply  uses  the  a  priori  belief  that  the  model  should  be  simple 
(smooth)  -  it  is  somewhat  ad  hoc  since  it  is  not  known  what  degree  of  derivative  is  appropriate  for  the 
model,  one  is  chosen  that  yields  a  visually  appealing  result  There  is  an  alternative  approach  that  can 
be  taken  when  sufficient  prior  information  exists  that  indicates  what  covariance  should  be  present  in 
the  model.  The  matrix  R^R  is  a  roughening  operator  which,  when  applied  to  the  model  vector  within 
the  functional  (equation  C.17)  yields  a  smoother,  or  more  correlated  model.  Tliere  exist  non-singular 
correlation  operators  that,  when  inverted  and  applied  to  the  model  vector  within  the  functional  C.16.  also 
yield  a  more  correlated  model.  We  denote  such  an  operator  by  If  we  return  to  C.I7  and  replace 

R^R  with  we  redefine  the  model  norm  by: 

llmi:*  =  (3.24) 


The  trade-off  functional.  U,  becomes: 

U  =  m'^C„lm  +  p  '{(Gin-d)’'C„,l(Gm-d)-X^}  (3.25) 

Following  the  steps  taken  above  to  define  equations  C.22  to  C.23  we  find  that  the  estimate  of  the  model 
can  be  cast  in  the  following  form: 

th  =  (cj„  ^  p  'CJCJ,C.)  '  n  ’G^C„^d  (3.26) 
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If  the  model  vector  is  due  to  a  stationary  process  and  if  Cmm  is  the  covariance  matrix  that 
describes  this  process,  this  operator,  when  applic'd  in  this  manner  (Equation  C.24),  will  introduce  the 
proper  degree  of  correlation  into  m.  Arguments  can  be  made  that  the  scattered  wavelields  should  be 
stationary  (ex.  Frankel  and  Wennerberg,  1987).  These  arguments  arc  ha.sed  on  the  assumption  that  the 
underlying  scatterers  are  distributed  randomly.  The  current  study  is  deterministic  in  nature,  the  number 
of  the  scatterers  is  not  known  beforehand  but  is  assumed  to  be  small.  The  scattered  waveficlds  we  are 
inverting  for  are  not  expected  to  be  stationary.  Since  wc  believe  the  model  we  seek  is  a  transient,  it  is 
probably  safest  to  consider  C,n,n  simply  as  a  smrxiihing  operator  that,  hopefully,  can  be  defined,  given 
a  priori  infonnation,  to  introrlucc  the  proper  degree  of  correlation  inlo  tlic  bansiciit  model  vector.  If 
we  follow  the  hypothesis  that  the  liKally  .scattered  phases  arc  generated  by  the  interaction  of  incident 
wavelields  and  It'pographic  teaturcs,  then  knowing  the  scale  length  of  the  totnigraphy,  and  the  velocity  of 
propagation  of  the  incident  energy,  wc  should  be  able  to  make  an  educated  guess  as  to  what  correlation 
should  be  present  in  the  model.  In  the  next  section  wc  will  asc  synthetics  in  an  attcmpit  to  determine  if 
this  Ls  true,  and,  if  so,  the  extent  to  which  the  result  can  be  improved. 

Tliis  result  (C.2.‘i)  is  somewhat  awkward  since  the  inverses  of  tlie  correlation  matrices  must  be 
computed  to  e,stiinate  the  model.  It  has  been  shown  (eg.  Tarantola,  1981)  that 


G'  ^G'C.  'CC . 


G'C„ 


(C. 


GC„,niG 


') = ( 


G^C,„lG^C,l)c,n,nG'  (3.27) 


given  th.it  C.,.;  -  C.C„u,S’^  and  G^  C,„5G  +  arc  positive  definite.  Since  C„,„,  and  have 
riK’plitc  lorrr;  itiey  are  positive  delinite  and,  thus,  so  arc  their  inverses.  .\s  a  rc.ojil. 


(g^CJG  .  C„.:„)g'c„^.  C . G''(C,...  ^GC„„„G')  '  (.3.28) 

and  equation  C.29  becomes; 

m  ..  C,„.,.G'  (/rC,.,.  .  GC.„„,G'  )  '  d  (3.29) 

II  tlie  model  and  noi.se  me  station.u'y,  stochastic.  priKcsses  ;ind  .ire  uncorrelatcd  with  each  other,  this 
result  has  the  form  of  the  stixhastic  inverse  of  Franklin  (1970),  It  is  only  nccessarv'  to  compute  a  single 
matrix  inverse  to  obtain  an  estimate  of  the  model.  However,  given  die  definition  of  G  in  Equation  C..^, 
this  m:itnx  is  larger  than  the  one  ajipiuinng  in  eqii.ilion  C.2.3  (N,(  by  N,i  vs  N,,,  by  N,„)  and  it  is  still 
nec.os.iry  to  invert  C„„  to  calcul.ite  mislit.  To  invert  for  nuxlcl  estimates,  equations  C.23  and  C.29 
lue  solved  in  steps  by  varying  //  ;uid  usirg  equtition  C.  19  to  determine  when  an  acceptable  misfit  has 
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hccn  achieved.  This  result  encompasses  the  Wiener  shaping  filter.  If  there  is  no  noise,  or  /i=0  (infinite 
emphasis  is  placed  on  misfit)  and  the  mexJel  is  assumed  to  he  uncorrelated  then  equation  C.29  simplifies 
to  equations  C.8  and  C.I4.  Hereafter,  we  refer  to  the  method  described  by  equation  C.29  as  methrxl  C. 


C.S.d  Wavelet  estimation 

In  Chapter  2,  we  found  that  the  best  practical  means  of  estimating  the  character  of  the  excitation 
function  (referred  to  as  the  wavelet  in  this  chapter)  is  via  beamforming.  This  estimation  procedure  is 
imperfect,  due  to  seismic  noise,  locally  scattered  phases  and  the  fact  that  the  incoming  wavefield  contains 
a  spectrum  of  wavenumbers.  This  beam,  which  will  be  denoted  by  b(i),  is  an  imperfect  representation 
of  the  wavelet  b(t),  but  is,  nonetheless,  assumed  by  the  preceding  theory  to  be  an  exact  representation. 
Tb  remedy  this  problem,  we  could  apply  the  approach  developed  by  Tarantola  (1982)  which  allows 
for  uncertainty  in  b(t).  Another  alternative  is  to  consider  the  current  problem  and  judge  whether  this 
additional  labor  is  likely  to  yield  much  improvement.  Assuming  that  the  signal  does  not  vary  between 
channels  but  the  noise  is  spatially  uncorrelated,  the  noise  in  6(f)  is  expected  to  be  a  factor  of  y/N  down 
from  the  noise  levels  in  the  individual  channels  (Lacoss,  197.5;  Husebye  et  al.,  1985).  The  NORESS 
array,  for  example,  has  25  sensors  so  the  noise  level  in  6(f)  should  be  roughly  ^  the  level  in  each  channel. 
In  the  next  section  we  attempt  to  assess  the  impact  of  such  an  error  on  a  synthetic  example. 

The  wavelet  deviates  further  from  the  uae  “excitation  function”  since  the  former  must  be 
truncated  at  a  reasonable  point.  The  impact  of  finite  wavelet  length  is  also  coasidered  in  the  next  section. 

C.4  Deconvolution  of  several  synthetic  time  series 

Prior  to  applying  the  technique  to  the  pre-processing  of  array  records,  the  method  and  software 
have  been  tested  on  a  variety  of  synthetic  datasets.  The  advantage  of  the  synthetic  data,  and  the  reason 
that  it  is  valuable  for  testing  purposes,  is  that  it  is  completely  controlled.  As  a  result,  all  aspects  of  the 
inversion  can  be  tested  since  it  is  known,  under  ail  circumstances,  what  the  outcome  should  be. 


C.4.a  Uncorrelated  noise  and  model 

The  initial  tests  have  been  made  on  a  simple  dataset,  one  in  which  the  wavelet  is  known  to  be 
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Figure  3.6  A  synthetic  minimuni  phase  wavelet  (described  in  the  text).  The  amplitude 
and  phase  specua  arc  displayed  in  the  lower  tigurc. 


a  simple  narrow-band,  time  decaying  sinusoid  dc.scribcd  by. 
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Figure  3.7  A  synthetic  time  series  produced  by  convolving  a  delta  function  (at  tM).6 
s)  with  the  wavelet  displayed  in  the  preceding  figure.  Uncorrelated.  Gaussian,  noise 
((t=0.(XX)2)  has  been  a^cd.  This  synthetic  h^  a  signal  to  noise  ratio  comparable  to 
the  best  events  in  the  NORESS  dataset. 


Th<  sinusoid  is  multiplied  by  the  hcavisidc  step  function.  H(t-to)  to  provide  an  abrupt  onset  at  t=to. 
n  '  >mtmls  the  rate  of  decay  of  the  exponential.  This  wavelet  is  the  result  of  an.  admittedly  simple- 
minded  attempt  to  mimic  a  typical  tcicseism  encountered  in  the  NORESS  dataset.  Although  the  input 
sinusoid  is  monochromatic  (at  frequency  fo).  as  illustrated  in  Figure  C.6  some  bandwidth  is  present, 
due  to  the  exponential  decay  and  Anile  length,  in  this  subsection,  all  input  time  series  are  constructed 
usn  g  uncorrelated  models  and  uncorrelated  noise  -  consequently  both  methods,  B  and  C,  give  exactly 
the  same  resulLs  (C„„  and  or  R)  .arc  diagonal  matrices.  The  time  scries  in  Figure  C.7  has 

re.s  ited  from  the  convolution  of  the  wavelet  with  a  delta  function  (at  t()=0.6  s)  and  the  addibon  of 
unc' irrelated  Gaussian  noi.se  (fT-=.0(X)2).  This  model  provides  an  interesting  challenge  to  methods  B 
anu  C  sitKe  these  approaches  both  employ  the  Lj  norm  to  assess  model  “size"  and  data  misAt  and  thus 
assume  thc.se  elements  obey  Gaussi...a  statisUcs.  In  this  synthcbc.  the  added  noise  belongs  to  a  stationary 
Gaussian  process,  the  model  parameters  clearly  do  not.  While  all  points,  except  one,  are  zero,  the  one 
exception  (the  d  funebon)  is  a  significant  outlier  (at  0.6  s  the  model  has  a  value  of  1.0).  The  output 
models  of  4  successive  inversions  arc  displayed  in  Figure  C.8:  the  data  misAts  are  displayed  in  Figure 
C.9.  From  inversions  I  to  3.  the  Lagrange  multiplier,  /i,  decreases  and  the  models  are  clearly  becoming 
rougher  since  more  emphasis  is  being  placed  on  misfit.  The  Lagrange  multiplier  acts  as  the  only  model 


Figure  3.8  Displayed  in  ihis  figure  are  mixlcls  resulting  from  deconvolving  ihe  wavelet 
displayed  in  Figure  C.6  from  (he  time  senc.s  displayed  in  Figure  C.7.  Fixir  inversions 
resulting  from  progressively  smaller  Lagrange  multipliers  are  included.  The  dashed  line 
in  the  upper  left-hand  figure  indicates  the  known,  correct  result.  The  lower  right-hand 
mixlel  has  been  selected  since  it  yields  an  ims  ikila  misfit  of  I.OI. 
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Figure  3.9  Di.splaycd  in  (his  figure  arc  data  residuals  resulting  from  deconvolving 
the  wavelet  displayed  in  Figure  C.6  from  the  time  scries  displayed  in  Figure  C.7.  The 
residuals  present  in  inversions,  resulting  from  the  four  models  displayed  in  Figure  C.8, 
arc  included.  The  lower  right-hand  figure  illustrates  thal  the  selected  model  provides 
overfitting  of  the  data  in  the  model  space  with  the  exception  of  the  points  in  the  vicinity 
of  the  delta  function  (at  0.6  s). 


(iV 


Figure  3.10  A  portkin  (rows  and  columns  from  150  lo  212  of  274)  of  Ihc  dala- 
rcsolution  matrix  (Rd)  for  the  inversion  displayed  in  figures  C.7  and  C.8.  The  model 
space  ends  at  row/column  175. 

srTKKithing  parameter  since  no  smoothing  is  caused  by  C„,m  and  (R^  R)  since  they  arc  identity  matrices. 
By  inspct-iion  of  equation  C.  17  it  is  clear  that  in  the  limit  where  //  vanishes,  we  get  the  model  with  the 
greatest  resolution.  From  inspection  of  Figure  C.8.  three  problems  wilh  the  preferred  model  (lower  right- 
h;iiul  comer)  immediately  become  apparent.  The  della  function  amplitude  is  underpredicted  (by  22‘^). 
Ihc  delta  function  is  accompanied  by  side  lobes  and  a  significant  amount  of  noise  has  been  iiitrcxluced 
into  the  mixlel.  In  Figure  C.9  it  can  be  seen  that  at  a  misfit  of  1.01.  the  most  significant  discrepancy 
between  input  and  predicted  data  occurs  at  the  location  of  the  delta  function  (0.6  seconds).  The  data  are 
being  over-fit  where  the  mtxlel  ha.s  support  (from  0.0  to  4.35  s)  and  Ihc  expected  amount  of  misfit  is 
(Kcurring  outside  tlie  mcxlel  sp;ite  (fnim  4.37  s  to  6.825  s). 

It  IS  possible  lo  describe  why  Ihc  data  misfits  arc  unusually  disUibuled  by  considering  the 
rn.inncr  in  which  the  predicted  ctila  d,,  relate  lo  Ihc  input  itila.  Using  equations  C.4  and  C.29  we  find 
Itial 
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Figure  3.11  The  upper  figure  illustrate.s  dependence  of  the  data  misfit  on  the  La¬ 
grange  multiplier  (fi)  for  the  inversion  considered  in  Figures  C.7  to  C.IO.  The  curves 
higl.'ighted  by  the  circles,  triangles  and  squares  represent  the  misfit  measured  over  the 
data  space,  the  model  space  and  the  portion  of  the  model  space  from  1.0  to  4.0  s.  The 
lower  figure  results  from  the  same  data  however  the  trade-off  illustrated  is  between 
data  misfit  and  model  si/e.  As  in  the  ujjpcr  figure.  //  decreases  from  left  to  right.  In 
both  ligurcs.  the  hori/.onlal  line  (at  an  rms  misfit  of  1.0)  represents  the  desired  misfit. 
The  input  model  lies  on  the  vertical  line  in  the  lower  figure.  In  the  lower  figure,  the 
horizontal  axis  scale  has  been  magnified  by  a  factor  of  14.4  to  allow  visual  inspection 
of  detail. 


Figure  3.12  The  left  ligure  illastralcs  the  model  that  yields  an  rms  data  misfit  of 
1.0  when  the  wavelet  (l  igure  C.6)  is  deeonvolved  from  a  time  series  identical  to  that 
displayed  in  Figure  C.7  except  that  the  noise  level  has  10  times  the  variance  (fT=0.002). 
The  higher  noise  level  has  resulted  in  a  smaller  main  lobe  with  larger  side  lobes.  The 
figure  on  the  right  illustrates  that  a  relatively  large  misfit  is  occurring  in  the  vicinity  of 
the  delta  function,  and  overfitting  is  occurring  over  the  rest  of  the  model  spaee. 


If  the  matrix  that  premultiplies  the  input  data  vector,  the  data  resolution  matrix.  Rj,  is  equal  to  I.  the 
identity  matrix,  the  data  are  predicted  exactly.  This  occurs  when  the  noise  has  zero  variance,  or  in  the 
lunit  where  the  Lagrange  multiplier  vanishes  (misfit  is  considered  to  be  of  the  upmost  importance).  In 
this  inversion,  neither  of  these  conditions  exist.  In  general,  GC„,„,G^  possesses  small  eigenvalues  and 
the  noise  correlation  matrix.  C„„,  is  needed  to  stabilize  tlie  inversion.  As  a  rcsull,  this  added  component 
prevents  the  data  re.solution  matrix  from  achieving  the  form  of  an  identity  mau^ix. 

The  impact  of  fiCn,,  is  greatest  at  the  lower  right  hand  comer  of  GC„,,„G^  where  values  are 
small  due  to  the  decay  of  the  wavelet.  T  tic  inverse  of  the  sum  /iGn„  -  GCn,,nG^  is  thus  significantly 
changed  from  the  inverse  of  GC,„,„G  ^  and.  as  a  result,  R,(  deviates  from  1.  In  this  case,  the  edge  effect 
appears  in  Figure  C.9  as  bloated  residuals  (data  misfits)  beyond  the  model  space.  A  portion  of  the  data 
resolution  matrix  for  this  inversion  is  shown  in  Figure  C.IO.  Within  the  mixlel  space  (at  row/columns 
less  that  176).  clearly  tlic  data  arc  well  resolved.  The  decay  in  power,  and  thus  the  loss  of  the  diagonal 
fomi,  begins  abruptly  at  row/column  176  -  the  first  point  outside  of  the  model  space.  The  periodicity 
in  ihe  energy  that  remains  beyond  this  limit  is  due  to  the  periodicity  in  b(t).  Examination  of  the  data 
resolution  matrix  reveals  that,  given  llie  rapid  decay  of  b(t),  it  is  not  possible  to  predict  stationary  data 
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Figure  3.13  The  upper  figure  illustralcs  a  synthetic  lime  scries  that  has  resulted  from 
the  convolution  of  1 5  delta  functions  with  lire  wavelet  introduced  in  Figure  C.6).  Un- 
correlatcd.  Gaussian,  noise  (rT=0.0002)  has  been  added.  The  magnitudes  and  locations 
of  the  delta  functions  are  revealed  by  the  solid  line  in  the  lower  figure.  The  dashed 
line  represents  the  mtxlcl  that  prodiKcs  an  rms  data  misfit  of  1.0. 


Figure  3.14  The  loll  figure  illustrates  the  model  that  yielils  an  mis  data  misfit  of  1.0 
when  the  wavelet  (Figure  C.6)  is  daonvolved  from  the  time  series  displayed  in  Figure 
C.7  after  adding  noise  to  the  wavelet  (5  the  amplitude  of  the  noise  in  the  data).  The 
misfit  is  illustrated  in  the  right-hand  figure.  This  figure  is  included  to  illustrate  what 
impact  error  in  the  estimation  of  the  wavelet  has  on  the  inversion. 


outside  (he  miKlel  space.  In  short,  the  noi.se  is  generated  as  a  stationary  process,  not  via  a  convolution, 
and  so  ihis  ds  cun  volution  prexedure  is  not  able  to  duplicate  it  -  there  is  an  edge  effect  made  noticeable 
bc-cause  of  the  rapid  decay  of  the  minimum  phase  wavelet.  A  gra[)hic  description  of  the  impact  of  the 
unevenly  distributed  data  residuals  on  the  chosen  model  is  given  by  trade-off  curves.  In  the  upper  half 
of  Figure  C.ll  is  displayed  the  trade  off  relationship  between  the  Lagrange  multiplier,  /z,  and  the  data 
mistil  The  three  curves  in  this  figure  were  computed  by  considering  mi.slit  over  the  full  data  space 
(circles),  the  mtxlel  space  (triangles)  and  the  {xirtion  of  die  model  space  from  1.0  to  4.0  s  (squares).  The 
solid  horizontal  line  in  both  figures  represents  the  desired  misfit.  Considering  either  misfit  over  the  model 
space,  or  the  data  space,  we  select  the  model  when  //  readies  a  valix'  of  since  at  this  point  the  misfit  lues 
reached  the  desired  level.  The  vertical  gap  between  the  misfit  assixiated  with  this  mixlel  and  the  misfit 
occurring  from  1.0  to  .1.0  s  (the  bulk  of  the  nuxicl)  illustrates  die  extent  to  which  ovcrlitting  is  tx'curring 
in  the  bulk  of  the  mtxlel  space.  Considering  only  the  misfit  occurring  in  the  model  space,  with  the 
exception  of  the  large  values  near  t=0.6  s.  the  nomiali/ed  misfit  is  roughly  0.29  -  a  significant  ovcrlitting 
IS  ixcurring.  To  a  large  cxleni,  die  overiiiting  that  ixcurs  in  most  of  the  model  space  compensates  for 
the  iinderliiting  that  (xcurs  near  the  delta  fuix'tion  at  O.h  s.  The  algorithm  has  trouble  fitting  the  delta 


Figure  3.16  A  syntholic  I'iik’  series  prcnlueed  convolving  a  delta  function  (at 
t=0.b  s)  with  the  wavelet  displayed  in  f-igure  C.6.  Correlated,  Gaussian,  noise  has  been 
prcxJuccd  by  convolving  unrorrclatcd  noise  with  a  Nixcar  function. 


function  and  thus  overfits  other  mcxlel  points  to  compensate.  In  liie  lower  figure,  we  display  the  trade-off 
between  data  misfit  and  model  size.  The  horizontal  axis  scale  has  been  magnified  by  a  factor  of  14.4  to 
allow  examination  of  detail.  As  a  result  of  this,  the  preferred  model  api^ears  further  from  the  bend  in  the 
trade-off  curves  than  it  actually  is.  The  sclccteii  mixkl  has  a  norm  of  roughly  0.65  -  the  input  model  is 
known  to  have  a  norm  of  l.G.  This  discrepancy  is  mainly  due  to  the  imderprcttiction  of  the  delta  function 
amplitude. 

A  second  test,  asing  the  same  input  model  but  adding  noise  w  iih  10  limes  the  variance  ((t*=.002) 
yields  the  output  model  and  data  misfits  pictured  in  Figum  C.12.  This  noise  level  is  roughly  equivalent 
to  the  worst  level  encountered  in  the  Scripps’  archive  ol  10.5  NORESS  tcle,seismic  recordings.  In  this 
noise  env  ironment,  at  a  misfit  of  !  01,  the  delta  function  is  still  located  in  the  cor  ec  place:  however,  its 
amplitude  is  judged  to  be  only  0.4  the  true  amplitude.  The  sidclobcs  have  increased  in  amplitude,  noise 
is  still  being  tit  into  the  mcxlel.  and  overfitting  is  occurring  only  in  the  model  space. 

In  Figure  C.I3  we  illustrate  the  degree  to  which  the  inversion  algorithm  is  aF'e  to  find  a  set  of 
15  delta  functions  convolved  with  the  minimum  phase  wave'.,.  (Figure  C.6)  when  the  result  is  added  to 
a  low  level  of  noise  (<7=0.0002)  The  input  data  arc  illustrated  in  the  upper  half  of  the  figure,  'he  result 
is  displayed  in  the  lower  half.  It  appears  ■‘lal,  given  a  !  iw  noi.se  environment,  the  approach  is  able  to 
uncover  a  l.irge  number  of  dclui  funefions  when  iliey  are  convolved  with  a  minimum  phase  wavelet.  To 
test  the  resist  tnee  of  '.he  inversion  .>chemcs  to  imperfect  estimation  of  the  wa.elcl,  b(l),  we  have  added 
uncorri'ated  Ci.-ussiaii  noise  (at  one-fifih  the  amplitude  present  in  the  data)  to  the  wavelet.  This  relative 
a.iplituuc  h.cs  been  chosen  to  si:  date  the  relative  error  expected  in  the  NORESS  data  set.  In  all  other 
respects,  we  have  duplicated  the  first  inversion  (/  c  ,  .  we  assume  the  wavelet  is  free  of  error).  Visual 
comparison  of  'he  result  tFigurc  C.14)  with  the  original  inversion  (Figures  C.7  .o  C.ll)  reveals  little 
degnuluion  of  model  quality. 

In  Figure  C.15  wl  truncate  the  w-vciet  after  20  points  (0.5  s).  As  expected,  large  “fruncation 
phases"  and  residual:,  have  resulfed  at  multiples  of  0.5  s  after  the  mixlel  power.  This  lest  rcprc.scnls 
a  severe  truncation  -  as  will  be  seen  laier,  recorded  lelcscismic  amplitudes  generally  decrease  rapidly 
enough  so  that  power  levels  at  the  end  of  the  window  are  low. 

The  problems  noted  above  may  be  due,  in  part,  to  the  adoption  of  the  /.2  norm  and  the 
underlying  assumption  of  Gaussian  sfatistics.  Mcihixls  B  and  C  arc  clearly  not  very  tolerant  of  (or 
resistant  to)  ctlicrs.  The  challenge  of  recovering  the  input  model  is  incrca.sed  when  the  model  and 
data  erri'i  possess  similar  covanances  (in  these  vkitasets  they  are  both  uncorrclated).  Fortunately,  it  is 
extremely  unlikely  that  tne  recorded  data  will  be  so  "long  tailed  "or  that  the  model  and  error  statistics 
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Figure  3.17  The  left  figure  illustrates  the  model  that  yields  an  rms  data  misfit  of 
1.0  when  the  wavelet  (Figure  C.6)  is  deconvolved  from  the  time  scries  displayed  in 
the  preceding  figure.  The  noise  was  assumed  to  be  uncorrelated  and,  as  a  result,  a 
significant  amount  of  the  correlated  noise  has  been  introduced  into  the  model. 


will  be  so  well  matched. 


C.4.b  Correlated  noise 

In  this  subsection,  we  consider  a  time  scries  produced  by  convolving  the  minimum  phase 
wavelet  with  a  single  delta  function  (at  0.6  s)  with  correlated  noise  added  to  the  result  (Figure  C.16). 
The  correlated  noise  was  generated  by  convolving  uncorrclatcd  Gaussian  noise  with  a  boxcar  function. 
In  Figure  C.  17  we  prc.scnt  the  preferred  model  that  resulted  from  using  methods  B  or  C  and  assuming 
the  noise  is  uncorrclatcd.  This  p<x)r  result  suffers  mainly  because  the  correlated  noise  has  been  fit  into 
the  model  -  the  algorithm  expected  the  noise  to  be  uncorrclatcd  and  thus  put  this  correlated  energy  in  the 
model.  In  the  next  figure  (C.18)  we  have  included  the  covariance  mauix  estimated  for  the  noise  using 
Equation  C.19  and  have  re-inverted  for  the  model.  Although  the  delta  function  has  acquired  significant 
side  lobes,  the  rest  of  the  model  is  relatively  free  of  noise  -  tlic  residuals  arc  well  distributed. 


C  4  c  Correlated  model 

In  this  subsection,  we  consider  a  more  physically  reasonable  model  -  a  correlated  transient 
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Figure  3.18  llic  left  tigurc  illustiates  the  niuclcl  that  yields  ;»n  rms  data  mislil  of  1.0 
when  the  wavelet  (Injure  C.h)  is  deeonMilved  from  the  time  scries  displayed  in  the 
preceding  figure.  The  true  noise  covariance  matrix  (C„„)  was  included  in  the  inversion 
causing  most  of  the  noise  to  be  excluded  from  the  model. 
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F'igiire  3.19  To  prixluce  the  suilfietu  time  series  (right  figure!  the  transient,  correlated 
nuxlel  (left  tieiirei  is  eomolved  \*i!h  thr  luinimiiin  phase  wavelet  displayed  in  Figure 
r /)  I'lKorr,  lated,  tiaii.-.'i.ui.  no.  e  i  t  ()otHl2)  h;-,,  tven  .idded. 


display..d  on  ifie  lelt  side  ot  Figure  t'  19.  lii  prixlixe  die  input  data  (nghi  side  of  Figure  C.  19)  uncorrclaied 
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Figure  3.20  The  left  figure  illustrates  (he  model  that  yields  an  rms  data  misfit  of  1.0 
when  the  wavelet  (Figure  C.6)  is  deconvolved  from  Jic  time  scries  displayed  in  the 
preceding  figure.  The  model  was  assumed  to  be  uncomclated.  As  a  result,  the  data 
misfit  (right  figure)  is  concentrated  at  the  location  of  the  model-power. 
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Figvu-e  3.21  The  left  figure  illustrates  the  model  that  yields  an  rms  data  misfit  of  1.0 
when  the  wavelet  (Figure  C.6)  is  deconvolved  from  the  time  series  displayed  in  Figure 
C.19.  The  Toeplitz  correlation  matrix  (Cm„),  estimated  by  applying  an  autocorrelation 
operator  to  the  known  model,  has  been  used  in  the  inversion.  The  data  residuals  (right 
figure)  are  well  (evenly)  distributed. 
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Figure  3.22  The  left  figure  illustrates  the  model  that  yields  an  ims  data  misfit  of  l.O 
when  the  wavelet  (Figure  C.6)  is  deconvolved  from  the  time  scries  displayed  in  Figure 
C.I9.  The  Toeplitz  correlation  matrix  (Cmm).  estimated  by  applying  an  autocorrelation 
operator  to  the  known  model  (and  approximating  the  result  witii  a  Gaussian),  has  been 
used  in  the  inversion.  The  data  residuals  (right  figure)  arc  well  (evenly)  distributed  with 
the  exception  of  a  slight  concentration  of  power  at  the  vicinity  of  the  model  power. 


noise  (a  =  0  00U2)  has  been  added  to  the  convolution  of  this  transient  with  the  minimum  phase  wavelet. 
In  Figures  C.2()  to  C.24  we  pre.sent  a  scries  of  inversions  involving  this  data  that  were  computed  to  test 
the  relative  worth  of  methods  B  and  C.  In  Figure  C.20  we  use  method  C,  the  stalistical  approach,  to 
invert  for  the  underlying  model;  however,  we  assume  that  the  model  is  uncorrelated.  The  model  has  not 
been  resolved  with  great  fidelity  -  it  possesses  some  uncorrelated  noise  due  to  overfilting  in  the  model 
space  itiat  compensates  for  the  undertitling  that  has  occurred  at  the  Icxation  of  the  model  power. 

In  Figure  C.21  we  attempt  the  same  inversion;  however,  prior  to  the  inversion  we  have  applied 
the  autocorrelation  operator  to  the  known  model  and  included  this  information,  within  the  cofTel.ition 
operator  pari  of  the  inversion.  Dc.spite  the  slight  concentration  of  power  at  roughly  0.6  s,  the 

data  residuals  arc  extremely  well  distributed  -  the  model  power  is  not  being  underfitted.  It  appears  that  the 
statistical  approach  works  well  in  spite  of  the  transient  nature  of  tlu;  model.  To  achieve  the  result  displayed 
in  Figure  C.22  we  approximate  the  true  covariance  matrix  with  a  Gaussian,  that  provides  roughly  the  same 
smixithing.  and  invert  to  lind  slightly  worse  residuals  -  the  transient  amplitudes  arc  underpredicted.  In 
Figures  C  23  and  C  24  we  use  method  B  and  define  the  roughening  matrix  to  be  discrete  first  and  second 
order  dilferenc  mg  operators  respectively.  Both  of  these  results  arc  clearly  inferior  to  those  achieved  using 
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Figure  3.23  The  left  figure  illustrates  the  model  that  yields  an  rms  data  misfit  of 
l.O  when  the  wavelet  (F.gure  C.6)  is  deconvolved  from  the  time  series  displayed  in 
Figure  C.19.  The  model  has  been  solved  for  using  equation  C.21,  a  first-derivative 
roughening  matrix  was  used.  The  data  residuals  possess  a  more  marked  concentration 
of  power  in  the  vicinity  of  the  transient  model  power  than  is  possessed  in  the  previous 
two  inversions.  As  a  result,  more  noise  has  been  fit  into  the  model  (note  the  “chatter" 
in  the  model  (left  figure)  and  the  overfuting  in  the  model  space  (right  figure). 

method  C  with  the  known  correlation  operator  Cmm.  but  only  slightly  worse  than  the  results  achieved 
when  Cmm  has  been  approximated  with  a  Gau.ssian.  In  both  of  the  inversions  achieved  using  method 

B,  the  transient  is  not  well  predicted,  compensatory  overfitting  is  occurring  in  the  model  space,  causing 
“chatter”  in  the  model. 

C. 4.d  Deconvolution  of  full  array  recordings  of  synthetic  and  recorded  events 

At  this  point,  we  rclum  to  the  images  presented  in  section  C.2  to  determine  the  extent  to  which 
we  can  use  methods  A  and  C  to  deconvolve  the  .seismograms  and  recover  the  radial  resolution  in  the 
images.  Recall,  from  section  C.2,  that  the  image  in  Figure  C..^  originated  by  convolving  a  Semipalatinsk 
teleseism  (b,(f),  pictured  in  Figure  C.2)  with  the  synthetics  used  to  generate  Figure  C.l.  In  the  same 
manner.  Figure  C.5  was  generated  from  C.l  using  the  tclcseismic  earthquake  beam  (b,(t)  pictured  in 
Figure  C.4). 

In  Figure  C.2.S  we  present  the  result  of  using  method  C  to  deconvolve  the  25  channels  underlying 
Figure  C..1.  We  have  added  a  small  level  of  uncorrelatcd  noise  to  each  channel  to  stabilize  the  inversion. 
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Figure  3.24  The  left  flgurc  illustrates  the  model  that  yields  an  rms  data  misfit  of 
1.0  when  the  wavelet  (Figure  C.6)  is  deconvolved  from  the  time  series  displayed  in 
Figure  C.19.  The  model  has  been  solved  for  using  equation  C.21,  a  second-derivative 
roughening  matrix  was  used.  As  in  the  previous  inversion,  a  significant  underfitting  of 
the  data  is  occurring  at  the  location  of  the  model  power.  As  a  result,  more  noise  has 
been  fit  into  the  model  (note  the  “chatter"  in  tlie  model  (left  figure)  and  the  overfitting 
in  the  model  space  (right  figure). 


We  use  with  the  exception  that  it  has  been  truncated  after  4(X)  points  (roughly  40  s)  and  tapered 
by  a  cosine  taper  from  points  361  to  4{X).  The  data  vector  (799  points)  has  been  tapered  by  a  cosine 
icr  from  (xiinis  401  to  799  to  reduce  edge  effects.  The  model  contains  400  points  and  spans  40  s. 
i'liis  combination  of  vector  lengths  was  chosen  to  reduce  edge  effects  -  as  seen  earlier,  any  truncation 
error  should  occur  N(,  points  beyond  the  nuxlel  power.  Using  the  parlance  introducca  in  section  C.2,  this 
approach  has  redaced  the  power  of  the  tail  by  roughly  I,*!  dB.  The  azimuthal  resolution  (width  of  the  main 
lobe  I  has  been  recovered  almost  exactly  and  the  image  side  lobes  have  been  recovered  lo  a  .similar  extern. 
The  deconvolved  trace,  station  aO,  is  displayed  below  the  image  and  bears  a  striking  resemblance  lo  the 
original  channel  dtsplaycd  in  Figure  C.  1.  Fortunately,  considering  that  this  inversion  required  roughly 
1 1  minutes  of  CPU  time  on  the  SDSC  Cray  Y-MP,  it  appears  to  have  been  a  great  success.  In  Figure 
C.26  wc  present  the  result  of  applying  method  A  to  deconvolve  the  same  25  time  series.  This  result  was 
achieved  after  several  (loughly  10)  iterations,  when*  parameters  (such  as  the  water  level)  were  chosen  to 
yield  a  model  estimate.  The  visual  ap(x:arancc  of  each  intermediate  mixlel  led  lo  the  selection  of  a  new 
(wc  hope  belter)  set  of  parameters,  and  .so  on.  This  result  (Figure  C.26)  is  inferior  to  the  preceding  one 
(Figure  C.2.*')  for  several  rca.s(His.  First,  the  ringing  has  led  to  a  slight  widening  of  the  main  lobe  (in  a 
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Figure  3.25  The  result  of  the  deconvolution  of  the  tcicseism  displayed  in  Figure  C.2 
from  the  seismograms  represented  by  the  trace  in  Figure  C.3.  Perfect  deconvolution 
would  result  in  this  image  being  identical  to  that  displayed  in  Figure  C.l.  Displayed  in 
the  lower  figure  is  the  deconvolved  channel  aO, 


radial  sense)  and  some  loss  of  a/imuthal  resolution.  The  surface  wave  packet  has  been  smoothed  leading 
to  a  further  reduction  of  :i/imullul  resolution.  The  character  of  the  image  “side  lobes"  has  not  been  well 
recovered  (sec  Figure  C.l).  The  "tail"  has  been  signilicantly  reduced  (by  roughly  15dB).  It  is  likely  that 
a  better  result  ctxild  be  achieved  after  several  more  visual,  and  subjective,  iterations;  however,  based  on 
this  disappointing  result,  wc  have  chosen  to  abtindon  this  approach. 

The  second  lest  image  (Figures  C.4  and  C.5)  pre.senLs  a  greater  challenge  for  method  C  since 
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Fig\irc  3.26  The  rc.su!(  of  (he  (JeeonvolutuMi  of  the  tcleseism  displayed  in  Figure  C.2 
frcwi  (he  seisniograins  R-presented  by  the  trace  in  Figure  C,.^  u.sing  (he  .spectral  division 
technique.  As  in  the  previous  figure,  perfect  'Icconvoiulion  would  result  in  this  image 
being  identical  to  that  displayed  in  Figure  C.l.  Displayed  in  the  lower  figure  is  the 
deconvolved  channel  aO. 


the  wavelet  is  more  protracted.  Using  method  C  we  can  reduce  the  tail  by  roughly  10  dB  (Figure  C.27). 
There  is  a  slight  edge  effect  at  the  beginning  of  the  deconvolved  Inrc.  The  side  lobe  character  has  been 
leii'vered  the  ;uiniutfi;il  resolution  h.is  not  K'cii  degraded 


C.3  Conclusions 
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Figure  3.27  The  result  of  the  deconvolution  of  the  tcleseism  displayed  in  Figure 
C.4  from  the  seismograms  represented  by  the  trace  in  Figure  C.5.  As  in  the  previous 
two  figures,  perfect  dcconvoluuon  would  result  in  this  image  being  identical  to  that 
displayed  in  Figure  C.l  Displayed  in  the  lower  figure  is  the  deconvolved  channel  aO. 

In  this  chapter,  three  approaches  to  the  seismic  deconvolution  problem  have  been  considered. 
The  first,  method  A,  is  a  straightforward  spectral  division  technique  that  uses  pre-whitening.  The  second, 
method  B.  is  akin  to  the  Occam’s  Ra^or  approach  of  Constable,  el  a/.(I987).  The  third  technique, 
developed  in  this  chapter,  is  essentially  the  stochastic  inversion  scheme  of  Franklin  (1970)  and  allows  the 
intrtxluction  of  a  priori  information  regarding  noi.se  and  model  covariances,  into  the  inversion  process. 
These  appraaches  have  been  applied,  with  vivying  degrees  of  success,  to  several  synthetic  data'^ts.  The 
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spectral  division  technique  has  been  found  to  be  inferior,  primarily  because  of  its  ad  hoc,  highly  subjective 
nature. 

It  has  become  clear  from  the  analysis  of  synthetic  data  that  the  measure  by  which  techniques 
can  best  be  Judged  is  the  manner  in  which  the  data  misfits  arc  distributed.  The  most  succesful  techniques 
yield  evenly  distributed  residuals.  Using  this  yardstick  to  Judge  success,  methods  B  and  C  are  equally 
capable  of  inverting  data  when  the  model  is  uncorrclatcd.  Application  of  these  techniques  to  models 
which  possess  significant  outliers  generally  yields  uns  itisfaclory  results.  This  is  most  likely  true  because 
these  techniques  use  the  L2  norm  to  assess  model  si/e  and  thus  assume  the  model  is  a  realization  of 
a  Gaussian  prtKess.  The  problem  seems  particularly  acute  whr  the  model  and  noise  possess  similar 
covariances  -  it  is  not  clear  where  the  energy  in  the  data  should  be  placed  -  much  noise  is  assumed 
to  be  part  of  the  model.  Given  accurate  a  priori  information,  it  appears  that  method  C,  the  statistical 
approach,  is  capable  of  the  best  results,  or  the  most  evenly  distributed  residuals,  when  the  input  model  is 
a  correlated  transient  pulse  -  the  most  physically  rca.sonablc  model  considering  that  we  are  most  interested 
in  deconvolving  seismic  records  for  surface  wave  packets  produced  by  scattcrers.  This  is  true  even  when 
a  Gaussian  Tocplitz  operator  is  substituted  for  the  model  covariance  matrix. 

Analysis  of  synthetic  images  created  by  applying  the  migrauon  technique  (Chapter  2)  to  data 
constructed  by  convolving  full  array  “recordings"  of  synthetic  scattered  waveticlds  with  real  telcseisms 
reveals  that  method  C  is  capable  of  significantly  enhancing  the  temporal  resolution  of  the  time  series  and 
thus  the  radial  resolution  of  the  images. 
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